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1  Introduction 

Over  a  number  of  years  HR  Wallingford  and  the  USACE  have  engaged  in  collaboration 
relating  to  levee  performance  and  wider  flood  risk  management.  In  particular  there  are 
significant  synergies  in  the  research  developments  in  both  the  US  and  UK  focused 
towards  the  improved  risk-based  management  of  levee  systems. 

In  September  2010,  USACE  commissioned  HR  Wallingford  to  undertake  a 
collaborative  project  to  explore  the  utility  of  the  tools  and  techniques  that  have  been 
developed  at  HR  Wallingford  and  widely  used  by  the  Environment  Agency,  within  the 
US  context.  Through  pilot  application,  the  aim  was  to  provide  a  forum  for  detailed 
scientific  exchange  between  lead  researchers  at  HR  Wallingford  and  in  the  USACE. 
This  report  presents  the  findings  of  this  study. 

The  primary  objective  of  the  project  was  to: 

•  Support  the  USACE  in  developing  and  implementing  robust  and  credible  risk- 
based  asset  management  analysis  tools  and  techniques  in  practice. 

Secondary  objectives  of  the  research  were  to: 

•  To  exchange  information  relating  to  asset  management  science  and  capabilities 
in  the  US  and  UK 

•  To  explore  the  applicability  of  HR  Wallingford  tools  and  techniques  in  the 
context  of  the  US  setting 

•  To  undertake  a  “gap  analysis”  to  identify  a  future  research  needs 

•  To  outline  a  programme  of  future  co-ordination  and  collaboration  with  USACE. 

The  project  was  organised  as  a  series  of  tasks: 

Task  1  Exchange  workshops  -  Providing  a  forum  to  exchange  asset  management 
science  and  capabilities  in  the  US  and  UK. 

Task  2  Case  study  and  application  -  Providing  a  pilot  application  of  the  HR 
Wallingford  tools  and  techniques  and  to  explore  their  utility  in  the  US  setting. 

Task  3  Gap  analysis  -  Providing  an  analysis  of  priority  research  needs  based  on  the 
case  study  and  on-going  research  elsewhere  within  ERDC/HEC. 

Task  4  Future  programme  of  research  and  development  -  Outline  a  forward 
programme  of  R&D. 

To  date,  two  exchange  workshops  have  been  held  (Task  1).  The  first  took  place  at  HR 
Wallingford  in  January  of  2011.  The  second  took  place  at  the  HEC  office  in  Davis, 
California  in  July  of  201 1.  Summary  information  regarding  the  workshops  is  provided 
in  Appendix  1.  The  primary  outcomes  of  these  workshops  were  the  agreement  of  the 
pilot  site  location  -  St.  Paul  Minnesota  and  agreement  on  the  models  that  were  of 
interest  -  FRE,  HR  BREACH  and  RELIABLE. 

This  report  describes  details  of  the  models  that  have  been  transferred  to  USACE  and 
their  application  on  the  pilot  site  (Task  2).  It  also  describes  the  details  of  a  gap  analysis 
(Task  3)  and  includes  recommendations  for  further  activities  (Task  4). 
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2  Summary  description  of  models 


2.1  RELIABLE 

Fragility  curves  that  describe  the  reliability  of  flood  defences  conditional  on  extreme 
hydraulic  loading  conditions  is  a  concept  that  originated  within  the  US,  USACE  (1996) 
and  has  been  adopted  by  the  UK.  These  curves  form  a  primary  input  to  risk  analysis 
models  like  FIEC  FDA,  EIEC  FRM  and  FRE  and  have  been  the  subject  of  recent 
research  both  in  the  US  and  UK,  Schultz  et  al.  (2010),  Simm  et  al.  (2009) 

The  Reliability  Tool  (RELIABLE)  is  a  prototype  software  tool  that  can  generate  levee 
or  flood  defence  specific  fragility  curves.  The  tool  was  developed  within  the 
FLOODsite  and  FRMRC  research  projects  (FLOODsite,  2007;  FLOODsite  2008).  It  is 
used  to  generate  structure-specific  fragility  curves,  based  upon  a  reliability  analysis  of 
multiple  potential  failure  modes.  The  logical  rules  that  link  the  failure  mechanisms  are 
represented  by  fault  trees.  The  ability  to  assess  the  reliability  of  a  defence  in  as  much 
detail  as  the  engineer  requires  is  an  important  feature  of  more  local  risk  analysis  and 
RELIABLE  has  the  potential  to  offer  this  capability.  It  is  currently  a  prototype  software 
tool  but  could  be  developed  into  a  robust  software  system  to  support  asset  management 
and  flood  risk  activities.  Further  information  on  RELIABLE  is  provided  in  Appendix  2. 


Figure  2.1  Example  fault  tree  and  RELIABLE  interface 
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2.2  HR  BREACH 

A  primary  aspect  of  flood  risk  analysis  relates  to  flood  defence  failure.  The  initiation  of 
failure  is  handled  probabilistically  through  the  fragility  curves.  Another  aspect  is  the 
development  of  breaches  and  the  subsequent  consequences  associated  with  inundation. 
Whilst  FRE  incorporates  a  conceptual  physics-based  approach  for  rapid  estimation  of 
breach  inflow  volume,  which  is  useful  for  broad-scale  modelling,  it  may  be  appropriate 
in  some  cases  to  undertake  a  more  rigorous  physics-based  breach  approach  locally  and 
utilise  the  outputs  within  FRE. 

Over  the  past  15  years,  HR  Wallingford  has  undertaken  development  of  a  breach  model, 
called  HR  BREACH  (See  Figure  2.2).  The  HR  BREACH  model  predicts  breach  growth 
through  flood  embankments  of  different  material  types  and  construction,  combining 
hydraulics,  soil  mechanics  and  structural  analysis  into  a  single  breach  prediction  model, 
Mohamed  et  al  (2002).  The  main  features  include: 

•  breach  development  through  piping  and/or  overflow; 

•  simulation  through  both  cohesive  and  non  cohesive  soils; 

•  simulation  of  homogeneous  or  zoned  structures,  also  including  grass  or  rock 
embankment  surface  protection; 

•  breach  growth  simulation  by  surface  soil  erosion  and  head  cut  with  discrete 
block  failure  processes  (including  the  choice  of  using  different  sediment 
relationships  for  different  embankment  types);  and 

•  probabilistic  approach  options  to  allow  the  user  to  simulate  material 
variability/uncertainty,  including  full  Monte  Carlo  simulation. 


Figure  2.2  GUI  of  HR  BREACH 


The  HR  BREACH  model  was  used  during  the  EC  IMPACT  (www.impact-proiect.net; 
2001  -  2004)  project  to  compare  predictions  against  field  and  laboratory  test  data.  The 
model  performed  best  (on  average)  in  comparison  with  a  number  of  other  breach 
models  from  across  Europe  and  the  US. 

The  European  FLOODsite  project  (www.floodsite.net)  ran  from  2004  -  2009  and 
offered  an  opportunity  to  further  develop  the  HR  BREACH  model.  The  research 
programme  allowed  for  a  more  rigorous  analysis  of  the  IMPACT  project  data,  and 
implementation  of  new  development  work  to  improve  model  performance.  Further 
information  on  HR  BREACH  is  provided  in  Appendix  3. 
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The  HR  BREACH  model  has  been  also  used  on  numerous  consultancy  and  research 
studies  by  HR  Wallingford  over  the  last  decade.  On  the  pilot  study  described  below,  the 
HR  BREACH  model  was  used  to  derive  improved  breach  inflow  estimates. 

2.3  FLOOD  RISK  ESTIMATOR  (FRE) 

FRE  is  a  decision  support  toolset  for  quantifying  economic  and  social  impacts  of 
flooding  for  present  day  conditions,  future  scenarios  and  different  flood  risk  mitigation 
measures.  It  integrates  information  on  fluvial  and  coastal  hazards,  with  information  on 
the  physical  flood  risk  system  (assets,  topography),  and  information  on  people  and 
property  in  the  floodplain,  to  quantify  the  overall  flood  risk. 


Figure  2.3  Interface  of  FRE 

FRE,  by  its  design,  is  flexible  and  therefore  capable  of  supporting  a  wide  range  of  flood 
risk  management  decisions.  Further  guidance  on  its  use  is  provided  in  the  supporting 
User  Manual,  HR  Wallingford  (20 10). 

FRE  incorporates  the  risk  analysis  approach  described  by  Gouldby  et  al  (2008).  The 
primary  components  are  extreme  value  distributions  to  describe  the  hydraulic  loads, 
fragility  curves  to  describe  the  reliability  of  the  system  flood  defences  and  depth 
damage  functions  to  evaluate  economic  consequences.  It  can  be  used  to  support  flood 
risk  management  support  decisions  relating  to,  for  example: 

■  What  is  the  existing  and  future  risk? 

■  Where  are  the  high  and  low  risk  areas? 

■  Which  assets  contribute  most  to  risk? 

■  What  are  the  key  drivers  of  risk  (probability  and  consequences)? 

■  Which  mitigations  measures  are  most  effective  at  reducing  risk? 

The  primary  outputs  of  FRE  are: 

■  Spatial  distribution  of  risk  expressed  as  Expected  Annual  Damage  (EAD) 
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■  Spatial  distribution  of  inundation  probability 

■  Levee  contribution  to  risk 

The  methods  used  in  FRE  have  been  applied  to  establish  the  National  Flood  Risk  in 
England  and  Wales  over  the  last  four  years  and  numerous  regional  and  local  studies. 
The  Environment  Agency  of  England  and  Wales  are  currently  being  trained  in  the 
application  of  this  approach  for  their  own  “in-house”  capability  in  running  the  model. 
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3  Application  to  St  Paul  pilot  site 


3.1  INTRODUCTION 

The  USACE  have  been  engaged  in  flood  risk  management  activities  within  the  St.  Paul 
area  of  Minnesota  for  many  years  and  they  are  currently  undertaking  analysis  on  this 
area.  This  area  was  therefore  selected  as  the  pilot  site  area.  The  following  description  of 
the  site  has  been  provided  by  USACE. 

St.  Paul  levee/floodwall  project  is  located  in  Ramsey  County,  MN  between  river  miles 
837.0  and  840.4  of  the  Mississippi  River.  The  original  project  was  authorized  by  the 
Flood  Control  Act  of  1958  and  construction  was  completed  in  1964.  Modification  to 
this  project  involving  a  raise  of  the  levee/floodwall  of  4  feet  was  authorized  by  the 
Water  Resources  Development  Act  of  1986  (PL  99-662).  Construction  of  this  portion  of 
the  project  was  completed  in  199x.  The  primary  purpose  of  the  project  is  flood  risk 
management.  It  includes  a  recreation  component  as  well  in  the  form  of  a  walking  / 
biking  trail. 

The  area  within  the  line  of  protection  is  approximately  445  acres  in  size  and  the  land 
use  is  primarily  commercial  and  industrial  in  nature.  It  is  located  entirely  within  the  St. 
Paul  city  limits  with  downtown  St.  Paul  across  the  river  to  the  north,  St.  Paul  aiiport  to 
the  east,  and  residential/commercial  neighborhoods  to  the  south  and  west. 

This  section  of  the  report  describes  the  application  of  RELIABLE,  HR  BREACH  and 
FRE  to  this  study  site. 

3.2  RELIABLE 

One  of  the  primary  inputs  for  the  flood  risk  model,  FRE,  is  the  fragility  curves 
associated  with  each  of  the  flood  defences.  Within  the  UK  a  generic  defence 
categorisation  system  has  been  developed  and  through  reliability  assessment  of  their 
different  failure  modes  a  generic  set  of  fragility  curves  has  been  developed, 
(Environment  Agency  (2007)).  These  curves  have  been  derived  for  broad  scale  analysis 
and  are  recognised  as  being  an  approximation.  This  approximation  is  necessary  for 
regional  and  national  scale  analysis  as  it  was  considered  impractical  to  undertake  a  site 
specific  analysis  for  all  defences  within  large  spatial  areas.  These  generic  curves 
comprise  approximately  60  different  classes.  Each  class  has  a  set  of  five  fragility 
curves  associated  with  different  condition  grades,  1  being  very  good  condition  and  5 
being  very  poor  condition. 

When  the  risk  analysis  is  being  undertaken  on  smaller  spatial  scales,  it  is  more  practical 
to  undertake  site  specific  reliability  analysis.  The  prototype  RELIABLE  model  has 
been  designed  to  undertake  site  specific  reliability  analysis  to  generate  reach  specific 
fragility  cures  for  use  in  flood  risk  analysis  models.  The  application  of  this  model  to  the 
St  Paul  Pilot  Site  is  described  below. 

During  the  PFMA  session  undertaken  by  the  USACE  and  stakeholders,  13  different 
modes  were  identified  as  potential  ways  that  the  St.  Paul  levee/floodwall  might  fail. 
Some  of  them  are  grouped  together  based  on  general  location  along  the  levee  alignment 
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and  are  represented  by  a  single  F1A  scenario  (six  failure  modes  are  represented  by  the 
“Middle”  F1A  scenario). 

Table  3.1  shows  a  list  of  the  identified  failure  modes,  extracted  from  USACE  2011  and  Figure 
3.1  shows  their  location.  The  defence  length  in  St  Paul  was  also  split  into  8  reaches  (see 
Figure  3.2). 

Table  3.1  Potential  failure  mechanisms  (from  USACE  2011) 


Performance  Failure  Modes 

PFM 

Description 

Station 

El 

Embankment  piping 

1+00 

E3 

Slope  instability 

135+00 

E5 

Embankment  overtopping 

140+00 

E6 

Piping  at  interface 

0+00 

PI 

Piping  associated  with  penetration  of  levee 

42+00 

P2 

Piping  associated  with  penetration  of  levee 

42+01 

W1 

1-wall  -  structural  failure 

85+00 

W3 

1-wall  -  global  stability 

85+00 

W5 

1-wall  overtopping 

85+00 

W12 

T-wall  overtopping 

65+00 

Cl 

Failure  to  install  closures 

160+00 

C4 

Failure  of  Fillmore  Ave  closure 

115+00 

11 

Failure  of  pump  station  (interior  flood  control) 

0+00 

PFMW1 
1 1  Wall  t 


I- WaB  Global  Stability 


T-wall  Overtopping 


PFMC4 

Failure  of  Earth  Closure  at  Fillmore 


*•  I +  M I ' '  2  y  . 

Piping  Associated  with  a  Penetration 


+>.  PFM1 1 

Failure  of  the  Pump  Station 


IJV  PFME5  OTP 
Overtopping  of  Embankment 


|  HPFMW5HI 
1-Wall  Overtopping 


PFMEI^PPPPPP 

Embankment  Piping  ^Station  0  to  2 


P||PFME6|P^ 

Piping  at  Interface 


EL  iJPVW  PFME3  ^ 

Slope  Instability :  Station  120  to  150 


‘  v 


h^Ppfmci  ru 

Failure  to  Install  Closures 


Figure  3.1  Potential  failure  mode  locations 
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Table  3.2  shows  the  details  of  each  defence  reach.  A  fragility  curve  was  provided  by  the 
USACE  for  each  defence  reach  (See  Figure  3.3). 


Figure  3.2  Levee  reach  locations 


Table  3  2:  Defences  data  in  the  St  Paul  study  area 


Reach 

Length 

(m) 

Elevation 

(m) 

Reach 

Type* 

Location 

1 

1,228 

219.34 

L 

Levee  from  Bluff  to  W.  Water  Street  Ramp 

2 

61 

219.71 

L 

Levee  from  W.  Water  Street  Ramp  to  Washaba 
Street 

3 

784 

218.75 

W 

1-wall  and  T-Wall  Riverfront 

4 

41 

218.27 

L 

Small  levee  section  behind  Comcast 

5 

81 

218.26 

W 

1  -wall  behind  Building 

6 

291 

218.25 

L 

Levee  with  closure  under  Lafayette  bridge 

7 

216 

218.13 

W 

Wall  with  pump  station  around  Fertilzer  Plant  with 
closure 

8 

2,234 

218.02 

L 

Levee  around  airport 

*L=Levee  -  W=Wall 
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St.  Paul  FPS 
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Figure  3.3  Reach  Fragilities  for  St  Paul  FPS  (from  USACE  2011) 
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To  assess  the  likelihood  of  defence  failure,  reaches  1,  3  and  8  have  been  selected  to  be  used 
within  RELIABLE.  Table  3.4  shows  the  assessed  failure  modes  for  each  reach. 


Table  3.4  Assessed  failure  modes  in  RELIABLE 


Reach 

Failure  modes 

Seepage  through  the  embankments  leading  to  piping  failure 

1 

Erosion  of  embankment  surface  by  overflow 

Sliding  of  the  embankment 

3 

Sliding  failure  of  the  wall  element 

8 

Erosion  of  embankment  surface  by  overflow 

A  number  of  modelling  runs  have  been  undertaken  with  RELIABLE  to  establish  the  fragility 
curves  for  reaches  1,  3  and  8.  Details  of  each  run  are  given  below. 

Reach  1 

For  Reach  1,  a  run  was  undertaken  to  obtain  the  fragility  curve  for  the  following  failure 
mechanisms  (See  Figure  3.3): 

■  Seepage  through  the  embankments  leading  to  piping  failure  which  covers  flow  through 
general  material  or  holes.  Predominantly  sandy  or  silty  fluvial  /  estuarial  dykes,  or  pre¬ 
damaged  clay  surface  layer.  Flows  /  pressures  lead  to  internal  erosion  and  piping 
(FLOODsite,  2007). 

■  Erosion  of  embankment  surface  by  overflow  which  can  lead  to  the  failure  of  the  crest 
and/or  the  rear  face  of  the  embankment.  If  the  flow  velocities  are  high,  grass  cover  (if 
exists)  can  also  be  damaged,  leading  to  direct  erosion  of  embankment  materials. 
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A  combined  run  for  the  above  two  failure  modes  has  been  undertaken  in  RELIABLE. 
Figure  3.4a  shows  a  comparison  between  the  fragility  curve  obtained  by  RELIABLE 
against  that  provided  by  the  USACE. 


Figure  3.4a  Reach  1  Fragility  curves  from  RELIABLE  and  USACE  2011 

For  completeness,  the  RELIABLE  and  USACE  2011  fragility  curves  were  compared 
with  those  that  would  be  used  to  represent  the  defence  reach  in  a  typical  UK  high  level 
risk  assessment  study.  Namely,  the  generalised  fragility  curves  associated  with  a 
Defence  Class  10  (Narrow,  turf  covered  levee  with  no  crest,  front  or  rear  protection). 
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The  generalised  curves  for  this  defence  class  are  shown  alongside  the  RELIABLE  and 
the  USACE  2011  data  in  Figure  3.4b 


Figure  3.4b  Reach  1  fragility  curves  from  RELIABLE,  USACE  2011  and  UK 
generalised  fragility  for  Defence  Class  10. 

From  Figure  3.4b  it  can  be  seen  that  the  generic  fragility  curve  for  condition  grade  3 
(fair  condition)  would  give  a  good  match  for  Reach  1.  It  should  be  noted  that  when 
undertaking  high  level  probabilistic  risk  assessments  in  the  UK,  if  a  defence’s  condition 
is  unknown,  grade  3  is  applied  with  greater  uncertainty.  As  well  as  giving  a  good 
approximation  of  the  fragility  curve  without  the  need  to  do  defence  specific  reliability 
analysis,  the  curves  could  be  used  to  explore  scenarios  of  investment  in  maintenance 
and  deterioration  by  allowing  the  risk  assessment  model  to  use  fragility  curves  which 
vary  through  the  simulation  time  from  condition  grade  1  (very  good)  through  to 
condition  grade  5  (very  poor). 

Reach  3 

For  Reach  3  two  runs  have  been  undertaken  to  obtain  the  fragility  curve  for  the 
following  two  failure  mechanisms  (See  Figure  3.5): 

■  Sliding  of  an  embankment  which  is  a  complicated  process  that  typically 
depends  on  an  adverse  combination  of  several  rather  than  a  single  factor.  The 
factors  that  could  lead  to  sliding  include  (FLOODsite,  2007): 

1-  Reduction  in  deadweight  of  organic  fill  material  in  embankment 
due  to  seasonal  desiccation. 

2-  Increase  in  hydrostatic  horizontal  loading  due  to  formation  of 
deep  tension  cracks. 

3-  Increase  in  uplift  pressures  beneath  embankment  in  confined 
founding  strata  that  is  in  hydraulic  continuity  with  flood  water. 

■  Sliding  failure  of  wall  element  which  is  a  simple  sliding  failure  of 
crown  wall  or  similar  element  driven  by  static  water  level  differences.  Most 
likely  after  other  effects  have  reduced  soil  strength.  Failure  occurs  when 
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horizontal  force  from  hydrostatic  water  level  difference  exceeds  net  shear 
strength  of  wall  /  foundation  junction  (FLOODsite,  2007). 

Figure  3.6  shows  a  comparison  between  the  fragility  curves  obtained  by  RELIABLE 
against  that  provided  by  the  USACE.  It  can  be  seen  that  the  USACE  2011  reliability 
assessment  shows  a  higher  likelihood  of  failure  given  load  than  the  RELIABLE 
assessment.  Or  put  another  way,  for  a  given  probability  of  failure,  the  load  is  typically 
one  foot  higher  in  the  RELIABLE  assessment  than  those  of  the  USACE  assessment. 
This  may  be  associated  with  additional  failure  mechanisms  that  may  have  been 
simulated  with  the  USACE,  2011  assessment  for  these  particular  1-Walls,  taking  on 
board  the  very  detailed  failure  mode  assessments  that  had  previously  been  undertaken. 
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Water  Level  (m) 


Figure  3.6a  Reach  3  fragility  curves  from  RELIABLE  and  USACE  2011 


The  RELIABLE  and  USACE  2011  fragility  curves  were  compared  with  the  UK 
generalised  fragility  curves  associated  with  a  Defence  Class  5  (Narrow,  brick  and 
masonry  or  concrete  wall  with  front  and  crest  protection).  The  results  are  shown  in 


Figure  3.6b 
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Figure  3.6b  Reach  3  fragility  curves  from  RELIABLE,  USACE  2011  and  UK 
generalised  fragility  for  Defence  Class  5 


It  can  be  seen  that  the  generalised  fragility  curve  for  ‘fair’  condition  (CG3) 
overestimates  the  performance  of  the  structure  when  compared  to  the  defence  specific 
assessments.  This  may  be  a  result  of  the  generalised  fragility  curves  being  unable  to 


EX6687 


Page  13 


R.  2.0 


Flood  Risk  Asset  Management 
USACE 


Wallingford 


represent  the  complex  nature  of  the  defence.  Of  note,  the  generalised  curve  for  poor 
condition  compares  reasonably  well  with  the  defence  specific  assessments. 

Reach  8 

For  Reach  8,  one  run  of  RELIABLE  has  been  undertaken  to  obtain  the  fragility  curve 
for  the  erosion  of  embankment  surface  by  overflow.  Figure  3.7a  shows  a  comparison 
between  the  fragility  curves  obtained  by  RELIABLE  against  that  provided  by  USACE 
2011.  Similar  to  Reach  1,  the  reliability  assessments  are  in  close  agreement. 

The  RELIABLE  and  USACE  2011  fragility  curves  were  compared  with  the  UK 
generalised  fragility  curves  associated  with  a  Defence  Class  10  (Narrow,  turf  covered 
levee  with  no  crest,  front  or  rear  protection)  and  similar  to  those  for  Reach  1 ,  a  very 
close  match  was  found  between  those  of  the  generalised  fragility  for  fair  condition  and 
those  of  the  more  detailed,  defence  specific  reliability  assessments.  The  results  are 
shown  in  Figure  3.7b 


Water  Level  (m) 


Figure  3.7a  Reach  8  fragility  curves  from  RELIABLE  and  USACE  2011 
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Figure  3.7b  Reach  8  fragility  curves  from  RELIABLE,  USACE  2011  and  UK 
generalised  fragility  for  Defence  Class  10 

The  following  points  may  be  concluded  from  the  reliability  modelling  work  undertaken: 

In  general  terms,  the  analysis  undertaken  with  RELIABLE  yielded  similar  results  to 
those  obtained  by  the  USACE.  It  is  of  note  that  these  results  were  also  in  broad 
agreement  with  those  that  would  be  obtained  from  using  the  UK  generic  fragility 
curves.  The  differences,  in  terms  of  risk,  that  arise  using  different  fragility  curves  are 
explored  further  below. 


3.3  HR  BREACH 

When  undertaking  flood  hazard  or  flood  risk  modelling  for  locations  that  are  protected 
by  flood  defences,  it  is  important  to  consider  the  presence  and  performance  of  the 
defences.  As  discussed  in  the  previous  text,  the  structural  reliability  of  the  defences  can 
be  represented  by  fragility  curves,  which  depict  the  probability  of  the  defence 
structurally  failing  when  subjected  to  a  given  load.  It  follows  that  defences  can  be  in 
one  of  two  system  states; 

•  Failed 

•  Non-failed 

In  the  non-failed  state,  the  defence  remains  structurally  intact  throughout  the  period  of 
raised  water  levels.  Any  flooding  that  may  occur  within  this  system  state  is  a 
consequence  of  either  water  levels  exceeding  the  defence  crest  level  leading  to  overflow 
conditions  or  exposure  to  waves  leading  to  overtopping.  Both  of  these  processes  can  be 
estimated  using  relatively  simple  relationships. 

In  the  failed  state  however  the  water  enters  the  floodplain  via  a  breach  in  the  defence 
which  may  begin  relatively  small  but  grow  over  time  due  to  the  inflowing  water  eroding 
the  defence  material  and  causing  instability  leading  to  collapse  of  the  exposed  faces. 
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With  the  breach  dimensions  varying  through  time,  one  needs  to  consider  more  than  the 
final  breach  size  to  estimate  the  inflow  accurately  and  the  flood  hydrograph  water  level; 
hence,  the  evolution  of  the  breach  during  the  flood  must  be  simulated. 

To  investigate  the  improvements  that  could  be  made  to  flood  estimation  in  the  St  Paul 
area  when  using  a  breach  simulation  rather  than  assuming  a  fixed  breach  condition, 
simulations  of  breaching  have  been  undertaken  using  the  numerical  model  HR  Breach, 
described  in  Section  2.2. 

The  breach  modelling  was  undertaken  for  Reach  1,  a  large  levee  constructed  of 
relatively  porous  sand  with  grass  protection.  The  principal  failure  mode  for  this  levee  is 
piping  leading  to  crest  lowering  and  subsequent  overflowing  and  erosion.  This  is 
followed  by  rapid  growth  of  the  breach  through  successive  phases  of  erosion  and  bulk 
failure  of  the  exposed  breach  faces.  Figure  3.8  shows  the  dimensions  of  the  levee  as 
represented  in  HR  Breach. 


29.25  ft 


Figure  3.8  HR  Breach  representation  of  the  Levee  cross-section  at  Reach  1 

The  modelled  flooding  conditions  simulated  were  for  a  1  in  1000  year  extreme  flow 
event  as  represented  by  Figure  3.9. 


EX6687 


Page  16 


R.  2.0 


Flood  Risk  Asset  Management 
USACE 


fcHR  Wallingford 


Stage,  Reach  1 


Figure  3.9  1  in  1000  year  stage  curve  for  Reach  1. 


The  first  few  simulations  of  breach  in  Reach  1  revealed  a  shortfall  in  the  HR  Breach 
functionality  which  prevented  it  from  accurately  representing  the  St  Paul  Floodplain; 
the  HR  Breach  model  previously  allowed  the  downstream  (landward  side)  flow 
conditions  to  be  described  by  either  a  cross-section  and  slope  (as  found  in  a  river  valley 
below  a  reservoir),  a  rating  curve  or  a  head  versus  time  curve.  The  downstream 
condition  influences  the  inflow  hydrograph  and  thus  the  breach  growth.  The  St  Paul 
floodplain  has  a  fixed  area  being  completely  bound  by  either  the  flood  defences  or 
rapidly  ground  along  the  southern  edge  of  the  floodplain. 

In  order  to  represent  the  St  Paul  assessment  better,  a  new  downstream  condition  was 
introduced  to  consider  the  finite  capacity  of  the  St  Paul  floodplain  basin  and  the 
associated  water  levels  on  the  floodplain  side  of  the  levee.  The  model  was  re-run  for  a 
number  of  different  scenarios  of  piping  leading  to  overtopping  and  erosion.  The  results 
are  presented  in  Figure  3.10.  The  scenarios  modelled  were  all  run  with  a  1  in  1000  year 
flood  hydrograph  and  are  described  as  follows; 

Scenario  1 :  A  one  foot  diameter  pipe  was  initiated  one  foot  above  the  toe  of  the  levee, 
as  the  water  level  reaches  the  pipe,  it  grows  until  instability  causes  the  pipe  to  collapse 
and  the  material  above  the  pipe  falls  plugging  the  pipe  but  lowering  the  crest  of  the 
defence.  The  levee  continues  to  protect  the  floodplain  until  the  lower  crest  is  exceeded 
and  overflow  commences,  eroding  the  levee  material  and  forming  a  breach.  The  breach 
grows  from  the  crest  to  the  bottom  of  the  levee  and  grows  in  width  until  the  water  levels 
in  the  river  fall  below  the  ground  level  at  the  levee. 

This  scenario  does  not  consider  the  depth  of  water  on  the  floodplain  side  of  the  levee 
when  calculating  the  evolution  of  the  breach  it  is  inappropriate  for  modelling  the  St  Paul 
site.  The  continued  growth  of  the  breach  gave  a  final  breach  width  of  306ft,  while  the 
maximum  inflow  discharge  was  20,600cfs. 
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Scenario  2:  The  same  parameters  and  chain  of  events  as  scenario  1,  however  this  time 
the  new  downstream  condition  is  used  to  simulate  the  finite  capacity  of  the  St  Paul  flood 
basin.  As  the  water  enters  the  floodplain  through  the  breach,  the  water  within  the 
floodplain  basin  rises.  This  reduces  the  water  level  differential  at  the  breach  and  thus 
reduces  the  inflow  rate.  The  new  downstream  condition  allows  much  better 
representation  of  the  breach  evolution  at  St  Paul.  The  inflow  rate  peaked  at  3,700cfs  and 
the  breach  width  was  36.4ft. 

Scenario  3:  The  same  parameters  and  chain  of  events  as  scenario  2,  however  this  time 
when  the  pipe  collapses,  all  of  the  material  above  the  pipe  is  washed  away,  commencing 
the  breach.  The  inflow  rate  peaked  at  2,200cfs  and  the  breach  width  was  37.3ft. 

Scenario  4:  The  same  parameters  and  chain  of  events  as  scenario  3,  however  this  time 
the  pipe  is  located  three  feet  below  the  crest  of  the  levee.  The  inflow  rate  peaked  at 
3,000cfs  and  the  breach  width  was  33.5ft. 

Scenario  5:  The  same  parameters  and  chain  of  events  as  scenario  3,  however  this  time 
the  pipe  is  located  half  way  up  the  levee  6.65ft  below  the  crest  of  the  levee.  The  inflow 
rate  peaked  at  2,300cfs  and  the  breach  width  was  32.9ft. 


Upper  Breach  Location  - 156700 


Time  (secs) 

Figure  3.10  Results  from  the  HR  Breach  simulations  of  a  breach  in  the  levee  at 
Reach  1  compared  with  those  from  the  USACE’s  HECRAS  breaching  simulation. 

The  analysis  has  shown  that  using  a  dedicated  breach  model  can  give  significant 
improvements  in  the  estimation  of  inflow  through  a  breaching  levee  than  is  achieved  via 
simple  assumptions  on  breach  width  and  height.  The  breach  model  simulates  the  growth 
of  the  breach  through  time  providing  better  estimates  of  breach  evolution  and 
dimensions  and  as  a  result  the  inflow  rate  evolves  with  the  breach.  The  model  allows 
scenarios  of  breaching  to  be  explored,  enabling  the  exploration  of  the  impacts  that  may 
be  expected  under  different  circumstances  of  breaching.  This  supports  the  use  of 
breaching  models  in  flood  risk  models  where  uncertainty  or  variability  can  be 
introduced  to  various  parameters  such  as  the  geotechnics  and  breach  initiation  for 
example.  In  summary,  the  breach  modelling  has  shown  that  significant  improvements 
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can  be  obtained  by  modelling  the  breach  rather  than  performing  the  traditional  approach 
of  considering  a  static  breach. 

3.4  FRE 

The  USACE  provided  a  range  of  information  that  was  necessary  to  set  up  the  FRE 
Model  for  St  Paul.  This  included: 

■  Digital  Terrain  Model  (DTM)  of  the  study  area 

■  Spatial  location  of  buildings 

■  Hydraulic  boundary  conditions 

■  Fragility  curves 

■  Spatial  location  of  levee's  and  gates  (closures) 

■  Depth-damage  information 

The  first  stage  of  the  analysis  involved  preparation  of  the  DTM  into  the  mesh  that  is 
used  by  the  inundation  model  (Figures  3.1 1  and  3.12). 


Figure  3.11  DTM  Ground  model  of  the  St  Paul,  Minnesota  Study  Area. 
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Figure  3.12  Computational  mesh  of  the  Inundation  model  for  the  study  area 


The  next  stage  involved  preparations  of  the  model  databases  relating  to  the  fragility  of 
each  levee  reach  and  the  spatial  location  of  the  properties  within  the  floodplain  area 
(Figure  3.13). 
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Figure  3.13  Spatial  location  floodplain  property 

To  determine  the  range  of  model  runs,  a  series  of  discussions  with  USACE  staff  were  held.  It 
was  noted  that  because  of  the  high  standard  of  the  defences  that  the  current  residual  risk  was 
relatively  low.  Therefore,  to  obtain  results  that  yielded  more  insight  into  the  performance  of  the 
underlying  models,  it  was  agreed  to  explore  the  outputs  with  different  assumed  crest  levels  that 
were  lower  than  reality.  Other  factors  that  arose  within  the  context  of  the  discussions  related  to 
the  difference  between  the  UK  generic  fragility  curves,  and  the  US  site  specific  ones, 
differences  between  the  UK  and  US  depth  damage  curves  and  the  resolution  of  the  model,  both 
in  terms  of  number  of  levee  reaches  and  DTM  resolution.  In  particular  it  was  noted  that  the 
pilot  site  was  relatively  small  scale  and  that  there  was  a  need  to  run  on  lager  spatial  scales  and 
hence  coarser  resolutions.  Hence  exploration  of  these  factors  was  of  importance.  With  regard 
to  the  number  of  levee  reaches,  further  insight  provided  by  USACE  identified  additional  reach 
lengths  that  could  be  defined;  sensitivity  analysis  was  therefore  undertaken  using  a  finer 
resolution  of  reach  lengths. 

A  summary  description  of  the  model  runs  that  have  been  undertaken  is  provided  in  Table  3.5.  It 
is  of  note  that  once  the  initial  model  has  been  set  up,  exploration  of  different  scenarios  is  not  an 
onerous  task,  generally  requiring  straightforward  modifications  to  one  or  more  of  the  database 
tables.  The  runtime  of  the  model  for  this  relatively  small  study  area  is  generally  less  than  1  min 
for  a  single  run. 
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Table  3.5  Summary  of  risk  model  runs 


Run  name 

DTM 

resolution 

(m) 

No.  of  levee 

reaches 

Fragility 

Depth 

damage 

functions 

Crest  levels 

Reference 

5 

8 

US 

US 

existing 

CL-1 

5 

8 

US 

US 

ex.  -lm 

CL-2 

5 

8 

US 

US 

ex.  -2m 

CL-3 

5 

8 

US 

US 

ex.  -3m 

UK  fragility 

5 

8 

UK  (generic) 

US 

existing 

UKdamage 

5 

8 

US 

UK 

existing 

Levee  resolution 

5 

54 

US 

US 

existing 

Gates_resolution 

5 

54 

US+UK 

(gates) 

US 

existing 

DTM  resolution 

50 

8 

US 

US 

existing 

A  summary  of  the  results  of  the  risk  model  are  provided  in  Table  3.6.  Selected  outputs  from  the 
model  runs  are  provided  in  Figures  3. 14  to  3. 1 7. 


Table  3.6  Summary  of  results  from  the  risk  model  runs 


Run  name 

Risk 

(EAD  £k) 

Risk 

(EAD  $k) 

%  change 
wrt  ref. 

Reference 

25.3 

15.9 

CL-1 

59 

37 

+132 

CL-2 

242 

152 

+856 

CL-3 

597 

374 

+2254 

UK  fragility 

26.0 

16.3 

+2 

UKdamage 

25.1 

15.7 

-1 

Levee  resolution 

25.5 

16.0 

+1 

Gates  resolution 

24.4 

15.3 

-4 

DTM  resolution 

30.1 

18.9 

+19 
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Figure  3.14  Spatial  distribution  of  residual  floodplain  risk,  EAD  ($k),  reference  case 
(existing  situation) 


Figure  3.15  Attribution  of  risk  to  levee  reaches,  EAD  ($k),  reference  case  (existing 
situation) 
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Figure  3.16  Attribution  of  risk  to  levee  reaches,  EAD  ($k),  global  crest  level 
reduction  of  3m  (CL-3) 


Figure  3.17  Attribution  of  risk  to  properties,  EAD  ($k),  global  crest  level  reduction 
of  3m  (CL-3) 
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The  FRE  risk  analysis  model  has  been  run  for  a  range  of  different  scenarios  to  explore 
the  sensitivity  of  the  model  to  different  inputs  and  to  demonstrate  the  range  of  potential 
outputs.  The  key  points  to  note  are: 

•  The  model  is  computationally  efficient  (less  than  a  minute  for  running  a 
single  scenario),  this  facilitates  exploration  of  a  range  of  different  “what 
if’  scenarios. 

•  The  model  considers  systems  of  flood  defences  and  the  potential 
consequence  of  flooding  associated  with  each  individual  element  in  the 
system. 

•  The  model  provides  a  range  of  outputs  that  can  support  risk 
management  decisions,  relating  to  maintenance  and  refurbishment  of 
existing  defences  or  strategic  planning. 

The  model  can  be  used  over  broad  spatial  scales  due  to  its  efficient  computational 
runtime.  This  runtime  is  achievable  as  a  number  of  simplifying  assumptions  are  made 
within  the  system.  These  assumptions  are: 

•  Volume  based  flood  spreading  algorithm  (ie  non  time-stepping, 
simplified  hydraulics). 

•  Independence  of  strength  of  different  levee  reaches. 

•  Dependence  of  hydraulic  load  within  a  flood  area 

•  Hydraulic  independence  of  separate  flood  areas 

Due  however,  to  the  simplifying  assumptions,  it  is  often  beneficial  to  deploy  the  model 
in  conjunction  with  more  detailed  hydraulic  analysis,  for  example.  It  can  also  be 
deployed  in  a  screening  role,  to  inform  further  modelling  studies  or  data  collation 
activities.  These  assumptions  do  however,  impose  certain  constraints  and  research 
initiatives  have  been  undertaken  and  are  currently  underway  to  address  these  issues. 
These  initiatives  are  discussed  further  in  Section  4. 
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4  Gap  analysis 


4.1  ONGOING  RESEARCH  INITIATIVES  OF  RELEVANCE 
Generic  fragility  curve  upgrade 

Because  of  the  heavy  use  of  generic  fragility  curves  in  broad  scale  flood  analysis  tools 
in  the  UK,  HR  Wallingford  is  currently  carrying  out  a  project  to  improve  these  curves 
for  use  in  various  asset  management  tools  and  to  fit  into  a  new  asset  classification  being 
adopted  by  the  Environment  Agency.  The  curves  will  take  account  of  latest  science, 
balancing  accuracy  against  ease  of  application  within  the  constraints  of  existing  flood 
risk  management  models.  The  project  also  links  to  doctoral  research  currently  underway 
on  sea  wall  fragility  curves. 

URBANFLOOD 

The  UrbanFlood  EC  FP7  project  (http://www.urbanflood.eu/Pages/default.aspx)  is 
creating  an  Early  Warning  System  framework  that  can  be  used  to  link  sensors  via  the 
Internet  to  predictive  models  and  emergency  warning  systems.  The  project  includes 
three  pilot  sites  to  apply  and  validate  at  full  scale  the  technology  being  developed  in  the 
project:  Amsterdam  (Netherlands),  Boston  (UK)  and  Rhine  River  (Germany).  The 
sensors  installed  at  the  various  sites  include  various  MEMS  modules  to  measure 
displacement  and  pore  pressure  and  fibre  optic  cables  able  to  detect  strains.  The 
gathered  data  are  used  for  dike  stability  evaluation  with  different  models  and  also, 
combined  with  an  Artificial  Intelligence  (Al)  component,  for  detection  of  anomalies  in 
dike  behaviour.  Detected  anomalies  trigger  assessment  of  the  likelihood  of  levee  breach 
and  the  consequences  in  terms  of  flood  propagation  and  damage  in  the  defended  urban 
area. 


AREBA  breach  model  for  system  risk  analysis 

The  HR  Breach  model  was  developed  around  a  decade  ago.  This  model  has  been 
widely  used  for  simulating  breaches  in  dams  and  embankments.  It  is  however,  too 
computationally  demanding  for  use  in  system  risk  analysis  models  that  are  applied  in 
practice.  A  simplified  breach  model,  AREBA,  has  therefore  been  developed.  Van 
Damme  et  al.  (2012)  (reproduced  in  Appendix  4,  for  ease  of  reference).  AREBA 
simulates  embankment  breach  processes  that  arise  as  a  result  of  erosion  from 
overflowing  water,  or  internal  erosion  due  to  pipes  that  are  formed  through  the 
embankment.  Discharge  through  the  breach  depends  on  the  breach  depth  and  width,  or 
pipe  dimensions.  AREBA  analyses  surface  erosion  failures,  headcut  erosion  failures  or 
piping  failures.  AREBA  has  been  validated  through  comparison  with  the  more  complex 
HR  BREACH  model  (Figure  4.1)  as  well  as  data  from  full  scale  breach  experiments. 
Full  details  of  the  model  and  its  validation  are  provided  by  Van  Damme  et  al.  (2012). 
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Figure  4.1  Comparison  of  the  outflow  discharge  from  a  levee  using  the  AREBA 
dynamic  breach  growth  model  and  the  well-established  HR  BREACH  Model 
Van  Damme  et  al (2012). 

Multivariate  extreme  value  analysis 

One  of  the  constraints  of  the  existing  FRE  methodology  is  the  handling  of  spatial 
dependence  within  the  flood  risk  analysis  model.  Within  the  UK  there  has  been  a 
significant  amount  of  research  into  statistical  methods  for  multivariate  extremes, 
Heffeman  and  Tawn  (2004)  and  this  has  been  investigated  for  use  within  the  context  of 
flooding,  Environment  Agency  (2011).  HR  Wallingford  has  developed  a  new  version 
of  FRE  that  builds  upon  this  research,  Wyncoll  and  Gouldby  (2012),  reproduced  in 
Appendix  5.  This  new  modelling  system  explicitly  captures  the  spatial  dependence  in 
extreme  floods  enabling  the  damage  associated  with  widespread  areas  in  a  single  flood 
event  to  be  evaluated  appropriately,  this  is  particularly  important  for  insurance 
purposes. 

RFSM  EDA  Inundation  model 

One  of  the  critical  aspects  of  system  flood  risk  models  is  the  simulation  of  inundation 
depths  and  velocities.  The  current  volume  based  approach  used  in  FRE,  is 
computationally  efficient  but  uses  simplified  hydraulics.  It  is  not  a  time  stepping 
dynamic  model  and  can  therefore  only  output  final  flood  depths.  In  risk  analysis 
models  the  inundation  component  must  be  computationally  efficient  because  of  the 
many  simulations  that  are  often  required  (100's,  if  not  1000's),  this  generally  exludes 
models  that  solve  the  full  shallow  water  equations.  HR  Wallingford  has  therefore 
developed  a  time  stepping  model,  RFSM  EDA,  that  is  computationally  efficient  yet 
sufficiently  accurate  for  flood  risk  analysis.  The  model  uses  a  relatively  coarse  grid  for 
flux  computations  but  utilises  the  detailed  resolution  topographic  information  through  a 
sub-grid  process,  which  ensures  it  maintains  its  accuracy.  The  model  has  been 
compared  against  a  range  of  other  inundation  models  where  it  has  been  shown  to 
perform  favourably  both  in  terms  of  accuracy  and  run  time,  Jamieson  et  al.  (2012). 

An  example  of  a  comparison  between  RFSM  EDA  and  a  full  Shallow  Water  Equation 
Model,  Infoworks  2D  is  shown  in  Figure  4.2  (also  see  Figure  4.3).  Table  4.1  shows  a 
comparison  of  runtimes  of  RFSM  EDA  against  other  models  on  a  series  of  benchmark 
tests  (Environment  Agency  (2010).  More  information  on  the  performance  of  RFSM 
EDA  is  provided  in  Appendix  5. 
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Figure  4.2  level  plot  for  RFSM-EDA  and  Infoworks  2D 
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Figure  4.3  Level/time  plots  for  other  inundation  models  for  the  same  test 
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Table  4.1  Comparison  of  runtimes,  model  timesteps  and  final  domain  volumes  (Jamieson 
et  al  (2012) 


Model 

Computation  time  in  minutes  for  each  test 

2 

4 

5 

8A 

RFSM-EDA  (mesh  A) 

0.015 

0.21 

0.23 

2.9 

Dynamic  RFSM 

0.19 

5.8 

9.8 

23.3 

TUFLOW  FV 

2.64 

24.5 

2.9 

72.6 

InfoWorks  ICM1 

0.73 

6.5 

0.7 

27.1 

JFLOW-GPU 

1.83 

2.3 

10.2 

16.2 

Lisflood-ACC2 

n/a 

1.97 

0.68 

n/a 

Fastest  other3 

0.4 

1.27 

0.6 

4 

Slowest  other3 

130 

282.8 

350 

307.8 

Note  1 :  The  runtimes  are  taken  from  InfoWorks  RS,  but  the  results  in  this  paper  are  from  InfoWorks  ICM; 
little  difference  is  expected.  Note  2:  Lisflood-ACC  runtimes  appear  in  (Neal  et  al.,  201  lb).  Note  3:  Fastest  and 
slowest  models  other  than  those  shown  in  this  paper,  but  appearing  in  Wright  et  al.  (2012), _ 


4.2  GAP  ANALYSIS  AND  PROJECT  RECCOMMENDATIONS 

Throughout  the  course  of  the  project  a  wide  range  of  topics  have  been  discussed  and 
potential  opportunities  for  future  collaboration  have  been  identified.  An  initial  list  is 
summarised  here.  This  list  is  however,  only  intended  to  initiate  discussion  at  the  final 
workshop  and  it  is  anticipated  this  will  be  amended  and  updated  following  the 
workshop 

Improved  RELIABLE  model 

Fragility  is  a  primary  requirement  for  the  FRE  model,  MDSF2  and  HEC  FRM  and  is 
embedded  within  the  Dutch  Reliability  model  HydraRing  and  FLOR1S.  The  existing 
RELIABLE  tool  is  prototype  only  and  further  investement  is  required  to  create  a  fully 
functional  and  robust  tool.  The  UK  Environment  Agency  have  indicated  some  interest 
in  further  development  of  the  tool  and  given  the  interest  in  these  approaches  from  the 
Netherlands,  there  maybe  an  opportunity  for  a  jointly  funded  project  to  develop  this 
tool. 

Multivariate  analysis  for  HEC  FRM 

The  HEC  FRM  model  requires  upstream  flows  from  different  river  systems  as  input.  It 
is  likely  that  flows  on  these  different  systems  will  often  be  partially,  not  fully, 
correlated.  The  Heffeman  and  Tawn  (2004),  method  can  be  used  to  analyse  flow 
records  from  different  gauging  stations  and  extrapolate  the  joint  probability  density  to 
extreme  values.  In  principle  this  can  be  applied  to  any  number  of  gauging  stations. 
Samples  that  form  the  input  to  the  HEC  FRM  model  can  then  be  generated  whilst 
preserving  the  appropriate  correlation  structure,  particularly  in  the  extremes.  HR 
Wallingford  staff  are  intimate  with  the  application  of  the  Heffeman  and  Tawn  (2004) 
methodology  in  the  context  of  flooding  and  would  welcome  the  opportunity  to  work  in 
partnership  with  counterparts  at  HEC  to  implement  this  system  in  HEC  FRM. 

Sub-grid  inundation  modelling  within  HEC  FRM 

Inundation  modelling  is  a  particularly  computationally  intensive  aspect  of  system  flood 
risk  models  and  the  computational  efficiency  of  the  inundation  model  can  be  the 
governing  factor  in  the  overall  perfomiance  of  the  system.  Both  HR  Wallingford  and 


EX6687 


Page  29 


R.  2.0 


Flood  Risk  Asset  Management 
USACE 


I  HR  Wallingford 


HEC  are  developing  sub-grid  inundation  models  to  speed  up  the  computation  whilst 
maintaining  accuracy.  There  may  well  be  merit  in  a  collaboration  to  share  ideas  and 
knowledge  with  the  overall  aim  being  to  implement  the  model  within  the  HEC  FRM 
system. 

Defence  system  state  importance  sampling  within  HEC  FRM 

The  total  number  of  simulations  required  of  the  inundation  model  governs  the  overall 
time  taken  to  complete  the  calculation  of  risk.  In  order  to  minimise  the  model  runtime 
it  is  therefore  desirable  to  limit  the  number  of  simulations  to  a  practical  number.  Both 
HEC  FRM  and  FRE,  currently  utilise  a  Monte-Carlo  sampling  method.  There  are 
however,  potential  opportunities  to  improve  the  computational  efficiency  of  both 
methods  by  implementing  Importance  sampling  (sometimes  known  as  variance 
reduction)  techniques.  Rather  than  sampling  at  random  from  the  range  of  defence 
system  states,  the  samples  are  stratified  and  weighted.  This  can  enable  a  significant 
reduction  in  the  number  of  samples  that  require  simulation.  Coupling  this  type  of 
sampling  method  with  an  efficient  sub-grid  computational  model  is  likely  to  yield  a 
significant  increase  in  the  computational  efficiency  of  HEC  FRM  and  FRE. 

Dynamic  Breach  growth  modelling 

The  rate  of  development  and  final  dimensions  of  breaches  in  levee's  can  have  a  critical 
influence  on  the  final  flood  extent  and  depths.  The  AREBA  model  is  a  computationally 
efficient  model  that  can  be  used  to  simulate  the  development  of  breaches.  Incorporating 
AREBA  within  HEC-FRM  will  enable  this  process  to  be  accurately  represented  without 
adding  additional  runtime. 
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5  Conclusions  and  Recommendations 

To  be  completed  post  workshop 
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Appendix  1  Workshop  Summaries 

Two  workshops  have  been  held  between  representatives  of  the  USACE  and  HR  Wallingford. 
The  first  in  January  2011,  HR  Wallingford,  the  second  in  July,  2011,  (HEC)  Davis,  California. 
The  workshop  agendas  and  outcomes  are  summarised  below. 
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HR  Wallingford  and  USACE  Workshop  : 

Flood  risk,  levee  safety  and  asset  management 

HR  Wallingford,  Tuesday  1st  February  -  Thursday  3rd  February,  2011 
Attendees 


USACE 
Noah  Vroman 
David  Schaaf 
Corby  Lewis 
Bob  Patev 
Neil  Schwanz 
William  Lehman 
David  Margo 
Alex  Roos 
Mitch  Laird 


HR  Wallingford 

Ben  Gouldby 
Caroline  McGahey 
Mark  Morris 
Mike  Panzeri 
Paul  Sayers  (S&P) 
Andy  Tagg 


Programme 


Workshop  Objective: 

The  objectives  of  the  workshop  are  summarised  under  two  components: 

1 .  Exchange  knowledge  on  levee  safety,  flood  risk  and  asset  management. 

2.  Develop  and  define  a  detailed  work-plan  (activities  and  programme)  for  implementing 
methods  and  software  tools  at  agreed  locations  to  produce  specific  agreed  outputs. 

These  activities  are  to  support  USACE  in  the  development  of  their  methodology  for 
prioritisation  of  maintenance  activities  for  risk  reduction. 

The  knowledge  exchange  aspects  will  comprise  topics  that  include,  but  are  not  limited  to: 

•  Condition  inspection 

•  Geotechnical  stability 

•  Closure  (active)  structures 

•  Fragility/reliability 
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•  Levee  deterioration 

•  Breaching  and  breach  modelling 

•  Risk,  uncertainty  and  sensitivity  analysis 

•  Economic  consequences  of  flooding 

•  Loss  of  life  from  flooding 

•  Hydraulic  modelling 

The  development  of  the  work-plan  will  seek  to  agree  specifics  relating  to: 

•  the  location  of  piloting  activities, 

•  the  software  tools  to  be  deployed  at  specific  locations 

•  data  gathering  activities 

•  site  location  visits 

•  further  workshop/s 

NB  -  The  programme  comprises  a  series  of  topics  for  discussion.  Presentations/slides  will  be 
provided,  it  is  however,  envisaged  these  will  be  informal  round  the  table  discussion  format  to 
stimulate  debate. 


Tuesday  lsf  February 

Morning 

9.00  Introductions  and  workshop  overview 

Overview  of  flood  risk  management  in  the  UK,  (HRW) 

Overview  of  flood  risk  in  the  US,  current  project  description  (USACE) 

Condition  inspection,  fragility  and  asset  management  in  the  UK  (HRW) 

Levee  monitoring  and  safety  in  the  US  (USACE) 


12:30-13:15  lunch 

13:15  -  14:15  Ship  Simulator  demonstration 
Afternoon 

14:15  (Description  of  Modelling  tools)  (HRW) 

System  risk  modelling 
Inundation  modelling 
Breach  modelling 
Consequence  modelling 
Reliability  modelling 
17:00  Finish 
Evening  meal  -  venue  the 


Wednesday  2nd  February 

Morning 

9:00  (Description  of  modelling  tools)  (USACE) 

Risk  analysis 
Hydraulic  analysis 
Consequence  analysis 
Reliability  modelling 
Levee  screening  approaches 
12:30-13:15  lunch 

13:15-14:15  Tour  of  the  physical  laboratory 
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Aftemoonl4:15(Case  studies) 

New  Orleans 
St  Paul’s  Minnesota 
Thames  Estuary 
Humber  Estuary 

17:00  Finish 
Evening  meal  -  venue  the 

Thursday  3rd  February 

Morning 

9:00  Software  demonstration 
Programme  development 
12:30-13:15  lunch 

Afternoon 

14: 15  (Summarising  programme  development  and  wrap  up) 

17:00  Finish 

Evening  meal  -  venue  the 


(USACE) 

(USACE) 

(HRW) 

(HRW) 


(HRW) 


(HRW/USACE) 
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HR  Wallingford  and  USACE  Workshop  : 

Flood  risk,  levee  safety  and  asset  management 
HR  Wallingford,  Tuesday  1st  February  -  Thursday  3ld  February,  201 1 


Minutes  (Final) 


Attendees: 


USACE 

Noah  Vroman  (NV) 

David  Schaaf  (DS) 

Corby  Lewis  (CL) 

Bob  Patev  (BP) 

Neil  Schwanz  (NS) 

William  Lehman  (WL) 

David  Margo  (DM) 

Alex  Roos  (AR) 

Workshop  Objective: 

The  objectives  of  the  workshop  are  summarised  under  two  components: 

•  Exchange  knowledge  on  levee  safety,  flood  risk  and  asset  management. 

•  Develop  and  define  a  detailed  work-plan  (activities  and  programme)  for 
implementing  methods  and  software  tools  at  agreed  locations  to  produce 
specific  agreed  outputs. 

These  activities  are  to  support  USACE  in  the  development  of  their  methodology  for 
prioritisation  of  maintenance  activities  for  risk  reduction.  It  is  currently  envisaged  the 
USACE  will  develop  a  framework  that  is  not  prescriptive  with  regard  to  modelling  tools 
and  implementation  but  provides  overarching  guidance  on  principles  and  concepts. 


HR  Wallingford 

Jonathan  Simm  (JS) 

Caroline  McGahey  (CM) 
Mark  Morris  (MM) 

Mike  Panzeri  (MP) 

Paul  Sayers  (PS) 

Andy  Tagg  (AT) 

Ben  Gouldby  (BG) 
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Summary  of  discussions  (Days  1  and  2) 

A  series  of  presentations  from  all  parties  took  place,  the  topics  were  wide  ranging  and 
included: 

•  Condition  inspection 

•  Deterioration 

•  Geotechnical  stability 

•  Closure  (active)  structures 

•  Fragility/reliability 

•  Levee  deterioration 

•  Breaching  and  breach  modelling 

•  Risk,  uncertainty  and  sensitivity  analysis 

•  Economic  consequences  of  flooding 

•  Loss  of  life  from  flooding 

•  Hydraulic  modelling 

•  Screening  tools 

•  Optimisation 

•  Life  cycle  modelling  and  continuous  simulation 

There  were  many  commonalities  between  the  approaches  used  and  in  particular 
common  definitions  of  risk,  use  of  fragility,  concepts  of  screening  tools,  condition 
inspection  and  uncertainty.  The  differences  lay  in  the  specific  practical  implementation 
and  these  were  largely  a  result  of  different  overarching  policy  and  regulation 
requirements. 

Specific  actions  arising  from  discussions  on  Day  1  and  2  are  summarised  below: 

Action  DM  to  provide  paper  on  levee  screening  approach  (and  manual?) 

Action  DS  to  provide  information  on  unified  toolbox  for  internal  erosion  and  piping 

Action  BG  to  provide  information  on  capping  of  damages 

Action  BG  to  provide  information  on  handling  of  deprivation  in  consequence  analysis, 
in  relation  to  WL’s  current  research. 

Action  HRW  to  circulate  draft  paper  on  continuous  simulation  to  WL. 

Summary  of  discussions  (Day  3) 

During  the  morning  of  day  3  software  demonstrations  took  place  in  particular  in  relation 
to  the  Wallingford  risk  analysis  and  decision  support  tool  and  the  Process  based  HR 
Breach  model. 

In  the  afternoon,  breakout  discussions  were  held  to  confirm  priorities  of  implementing 
tools,  details  of  programme  of  work  and  points  of  contact.  The  outcomes  of  the 
discussions  are  summarised  below. 

It  was  agreed  that  St  Paul’s  would  be  the  focus  for  the  pilot  site  analysis  and  trials.  It 
was  noted  however,  that  a  single  site  wouldn’t  cover  a  full  range  and  provide  a  full  test 
for  the  models.  An  option  to  overcome  this  issue  was  to  experiment  with  different 
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modifications  to  the  St.  Paul’s  site  (eg,  modify  defence  crest  levels  to  facilitate  more 
damage).  The  primary  model  applications  and  outcomes  of  interest  are  described  below: 

Risk  analysis  model: 

As  a  general  principle  it  was  agreed  that  the  main  focus  of  effort  should  be  on 
implementing  the  software  as  is  rather  than  make  extensive  changes  to  the  Wallingford 
pre-processing  tools  for  example,  to  be  compatible  with  USACE  data  sets. 

It  was  agreed  to  run  the  UK  risk  analysis  model  on  the  St  Paul’s  area.  The  standard 
outputs  the  HR  Wallingford  risk  analysis  model  currently  produces  are: 

•  Maps  showing  EAD  over  the  floodplain  area 

•  Maps  showing  floodplain  EAD  attributed  to  levee  sections 

•  Maps  showing  annual  probability  of  exceedence  over  the  floodplain  area 
During  the  discussions  it  became  apparent  that  loss  of  life  is  a  key  metric  for  USACE 
applications.  There  was  some  discussion  on  whether  to  include  HRW’s  Dynamic  (as 
opposed  to  static)  Rapid  Flood  Spreading  Model  (RFSM)  in  the  risk  analysis  model.  It 
was  felt  however,  that  this  may  compromise  runtimes  and  that  a  simple  topographically 
based  approach  would  suffice  for  the  current  study.  This  will  mean  additional  output 
will  become  available: 

•  Expected  life  loss  (ie  risk  of  life  loss)  over  the  floodplain  area 

•  Expected  life  loss  attributed  to  levee  sections 

Action  EIRW,  to  implement  simple  loss  of  life  approach.  WL  can  provide  input  and 
advise  on  USACE  LIFESIM  approaches. 

To  enable  comparison  with  HEC  FRM  approaches,  it  was  agreed  it  would  be  useful  to 
undertake  single  hydrodynamic  simulations  for  specific  realisations  (ie  single 
hydrographs  and  system  state  combinations).  Model  cascade  could  include  HEC-RAS, 
HR  BREACH  and  Dynamic  RFSM.  Action  HRW  to  consider  model  coupling  and 
advise  DM. 

It  was  felt  useful  to  undertake  the  risk  analysis  using  the  UK  “generic”  fragility  curves 
and  compare  the  results  with  the  site  specific  fragility  curves  that  have  already  been 
derived  by  USACE.  This  will  provide  information  that  is  useful  for  understanding 
issues  of  transferability  to  other  sites  (ie,  do  the  generic  curves  act  as  a  reasonable  first 
pass,  without  having  to  undertake  detailed  reliability  analysis  for  each  defence). 

Action  HRW  to  set  up  risk  analysis  model  on  the  St.  Paul’s  study  site  and  undertake 
analysis  as  described  above. 

Breach  process  model: 

It  was  felt  that  it  would  be  of  benefit  to  explore  the  application  of  the  Breach  process 
model.  Action  HRW  to  configure  breach  model  for  the  St  Paul  study  site  ready  for 
transfer  to  USACE. 

RELIABLE:  It  was  felt  that  it  would  be  of  benefit  to  apply  the  RELIABLE  model  to 
compare  with  fragility  curves  that  have  already  pr.  Action  HRW  to  configure 
RELIABLE  model  for  the  St  Paul  study  site. 

RAFT:  The  USACE  screening  tool  is  similar  to  the  UK  RAFT  model.  It  was  felt  there 
would  be  merit  in  comparing  the  two  tools  in  the  context  of  the  St  Paul’s  area.  Nb: 
HRW  has  developed  the  RAFT  tool  for  the  Environment  Agency  (EA)  and  its 
application  on  St  Paul’s  will  require  consent  from  the  EA.  Action  HRW  to  enquire 
with  the  EA  about  use  of  the  tool  on  St  Paul’s. 
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Full  uncertainty  analysis:  Whilst  this  was  an  area  of  interest,  this  was  considered  as  a 
lower  priority  than  tasks  above  and  will  not  be  applied  as  a  specific  priority  at  St. 
Paul’s. 

Automated  intervention  optimisation:  This  current  HRW  research  initiative  was  of 
interest,  but  it  was  felt  this  was  at  a  too  early  stage  and  of  less  relevance  for  the  Corp’s 
current  project.  This  will  not  be  implemented  at  St.  Paul’s. 


Timescales 

A  number  of  key  milestones  were  identified  together  with  associated  timescales: 

1 .  HRW  visit  to  St  Paul’s  -  The  primary  objectives  of  this  trip  are: 

a.  HRW  site  familiarisation 

b.  Presentation  of  initial  results 

c.  Transfer  and  training  of  software 

Date:  Provisionally  agreed  as  early  May  201 1  -  Action  BG  to  contact  DM  and  confirm 
dates. 

2.  USACE  visit  to  HRW-  The  primary  objectives  of  this  trip  are: 

a.  Familiarisation  of  approaches  for  wider  group  of  USACE  personnel. 

b.  Presentation  and  discussion  of  updated  results. 

c.  Troubleshooting  software 

d.  Final  software  transfer 

Date:  Provisionally  agreed  as  end  June  201 1  -  Action:  BG  to  contact  DM  and  confirm 
dates. 


3  Final  USACE  reporting:  August  2011. 


Communication 

It  was  felt  most  efficient  for  points  of  contact  (POC)  to  be  identified  and  direct 
communications  on  specific  technical  issues  to  occur,  keeping  BG  and  DM  cc’d. 
POC’s  are  detailed  below: 


USACE 

Fragility/reliability/breaching 
Hydraulics  and  hydrology 
Consequences 
Data 


David  Schaaf 
Corby  Lewis 
Will  Lehman 
Andrew  Sander 


HRW 

Fragility 

Breaching 

Risk  software 

Data/consequences 


Jonathan  Simni 
Mark  Morris 
Caroline  McGahey 
Mike  Panzeri 
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HR  Wallingford  and  USACE  In  Progress  Review  Meeting: 

St.  Paul  Pilot  Study 

Flood  risk,  levee  safety  and  asset  management 
Hosted  by  USACE,  Tue  26th  July  -  Fri  29th  July  2011 

Programme  (Draft) 

Meeting  Objective: 

The  objectives  of  the  meeting  are  summarised  under  three  components: 

3.  Review  application  of  risk  assessment  models  and  tools  for  the  St.  Paul  Pilot  Study 

4.  Demonstration  and  handoff  of  the  HRW  software  tools  to  USACE  for  the  St.  Paul  Pilot 
Study 

5.  Develop  and  define  a  work-plan  (activities  and  programme)  for  completion  of  the  pilot 
study  including  identifying  opportunities  for  further  collaboration 

These  activities  are  to  support  USACE  in  the  development  of  their  methodology  for 
prioritisation  of  activities  for  risk  reduction. 

The  review  of  risk  assessment  models  will  comprise  topics  that  include,  but  are  not  limited  to: 

•  Data  and  inputs 

•  Methodology  and  computations 

•  Outputs  and  conclusions 

•  Comparison  of  methodologies  (USACE,  HRW,  WAT/FRM) 

The  demonstration  and  handoff  of  software  tools  will  include 

•  Hands  on  use  of  the  tools  by  USACE 

•  Handoff  of  the  software  tools  and  input/output  files 

The  development  of  the  work-plan  will  seek  to  agree  specifics  relating  to: 

•  List  of  action  items  required  to  complete  the  pilot  study 

•  Further  meetings 

The  programme  comprises  a  series  of  topics  for  discussion.  Presentations/slides  will  be 
provided,  it  is  however,  envisaged  these  will  be  informal  round  the  table  discussion  format  to 
stimulate  debate. 
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HR  Wallingford 

Working  with  water 


Tuesday,  26th  July 

9.15  Depart  Hallmark  Inn  for  HEC  Office 

9.30  Arrival  at  HEC  Office  (Shewbridge) 

Overview  and  Introductions 

9.35  Overview  of  Sacrament  River  Flood  Risk  Management  System  (Tibbits) 

10.30  Overview  of  RD 1000  Levee  System  (Tibbits) 
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11.30  Q&A 

12.00  Depart  HEC  Office  for  RD1000  Levee  System 
Stop  for  Lunch  En  Route 

13.30  Tour  of  RD1000  Levee  System 

17.00  Linish 

Depart  for  Hallmark  Inn 

Happy  Hour  at  Hallmark  Inn 
Dinner (venue  tbd) 


Wednesday,  27th  July 

08.30  Data,  Model  Inputs,  and  Risk  Assessment  Lramework 

09.00  Discussion 

09.30  Preliminary  Model  Results 

10.00  Discussion 

10.30  Break 

10.45  Reliability  Modeling 

Condition  Grade  Methodology 
Reliable  Methodology 

11.30  Discussion 

12.30  Lunch 

13.30  Breach  Modeling 
HR  Breach 

14.00  Discussion 

15.00  Break 

15.15  Model  Variations  and  Sensitivities 
Overtopping  Risks  and  Breach  Risks 
Levee  Crest  Elevation 
Model  Resolution 

USACE  Lragility  and  HRW  Lragility 
16.00  Discussion 
17.00  Linish 

Happy  Hour  at  Hallmark  Inn 


(HRW) 

(HRW) 

(HRW) 

(HRW) 

(HRW) 
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Dinner (venue  tbd) 

Thursday,  28th  July 

08.30  Software  Demonstration  (HRW) 

09.30  USACE  Trial  Use  of  Software 
12.00  Lunch 

13.00  Overview  of  USACE  Risk  Assessment  for  St.  Paul  (Hauck/Schwanz) 

15.00  Break 

15.15  Comparison  of  USACE  and  HRW  methods  and  tools 
17.00  Finish 

Depart  for  Sacramento  for  tour 
18.00  Start  tour 

19.30  Diner  (venue  tbd) 

Friday,  29th  July 

08.30  HEC  Mission  and  Current  Activities  (Harris) 

09.30  HEC  WAT/FRM  Application  for  St.  Paul  (Baker) 

10.30  Break 
10.45  Discussion 
12.00  Lunch 

13.00  Develop  Work  Plan 
15.00  Finish 


List  of  HRW  Attendees 

Gouldby,  Ben 
McGahey,  Caroline 
Panzeri,  Michael 

List  of  USACE  Attendees 

Baker,  Penni 
Empson,  Bill 
Harris,  Jeff 
Hauck,  Kari 
Lehman,  Will 
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Lewis,  Corby 
Margo,  David 
Needham,  Jason 
Patev,  Bob 
Roos,  Alex 
Sander,  Andrew 
Schaaf,  David 
Schwanz,  Neil 
Shewbridge,  Scott 
Terry,  Tom 
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1  Introduction 

The  RELIABLE  model  has  been  developed  on  the  EU  funded  FLOODsite  Project.  The  prototype 
software  is  captured  in  four  main  components: 

■  Structure  specific  fault  trees  -  constructed  externally  by  the  user  within  OpenFTA  software. 

■  Limit  State  Equations  (LSEs)  -  comprised  within  one  or  more  Dynamic  Link  Library  (DLL) 
constructed  from  Fortran  subroutines,  and  based  upon  the  Task  4  failure  mode  report.  The 
software  does  however,  provide  the  opportunity  to  extend  the  library  of  LSE’s  if  this  is  a 
requirement. 

■  Uncertainties  on  the  input  parameters  of  the  LSEs-  input  through  a  spreadsheet  interface. 

■  Numerical  integration  -  Monte-Carlo  simulation  implemented  through  c#  code  behind  a 
Microsoft  Excel  spreadsheet  interface. 

These  components  are  captured  in  Figure  1.1.  The  primary  outputs  of  the  software  are  a  fragility 
curve  for  a  specific  structure  or  the  annual  probability  of  failure. 


Figure  1.1  Components  of  the  Reliability  calculator  software 

1.1  System  requirements 

To  operate  effectively,  the  reliability  calculator  requires: 

■  Windows  XP  or  later 

■  Microsoft  Excel  2003  (earlier  versions  may  be  acceptable  but  have  not  been  tested). 

■  Microsoft  .NET  Framework  2.0  (may  be  downloaded  free  of  charge  from 
http://www.microsoft.com/downloads/details.  aspx?FamilyID=0856eacb-4362-4b0d-8edd- 

aab  1 5  c5  e04f5 &DisplayLang=en 

■  Local  Administrator  rights  for  one  or  two  steps  of  the  installation  process. 

*  A  folder  C:\TEMP  with  write  access. 
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2.  Description  of  the  system  components 


2. 1  Limit  state  equations 

The  calculator  uses  a  DLL  that  is  constructed  from  a  library  of  LSEs  (each  LSE  is  a  separate 
FORTRAN  subroutine  (Figure  2.).  The  calculator  comes  provided  with  a  DLL  (Task7_LSEs.dll)  built 
from  over  60  LSEs  based  on  the  FLOODsite  Task  4  report  (Allsop  et  al,  2007).  An  example 
FLOODsite  Task  4  template  is  shown  in  below. 


Bal.5aiii  Uplifting  of  impermeable  layers  behind  earth  embankment 


Summary:  Uplifting  behind  embankments  occurs  if  the  difference  between  the  local  water  level  h,  and 
the  water  level  “inside”,  hb  is  larger  than  the  critical  water  level  hc 


Reliability  equation: 

The  reliability  function  is  expressed  by: 


2  =  m0  ’  hc  -  mh  •  Ah 

where: 

hc  =  critical  water  level  [m] 

Ah  =  difference  between  local  water  depth  in  front  of  dike  and  water  level  in  the  floodplain  [m] 

m0  =  model  uncertainty  factor  [-] 

mh  =  model  uncertainty  factor  for  damping[-] 


Loading  equations: 

Resistance  (strength)  equations: 

Ah  =  h-hb 

h  _  Y  wet  Y  w  ^ 

Yw 

Parameter  definitions: 

ywet  =  saturated  volumetric  weight  of  the  impermeable  soil  layers 
yw  =  volumetric  weight  of  the  water 
d  =  thickness  of  the  impermeable  layers 

h  =  water  level  on  the  river  [m] 

hb  =  water  level  in  the  floodplain  [m] 


Sources  of  failure  mechanism  equations  /  methods: 
Vrouwenvelder  et  al.  (200 1 ) _ 

Sources  of  uncertainties  in  failure  equations  /  input  parameters: 

Vrouwenvelder  et  al.  (200 1 ) 

Remarks: 
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REAL  FUNCTION  LSEBA1  5  A I  I  I ( VA L UE 5 , V AL UE 5  COUNT ,  E RROR NO,  E R ROR MS G) 

USE  LSESupport 

! Standard  Arguments  for  an  LSE  Function 

I  NTEGER,  I  NT E NT ( I  N )  :  :  VALUESCOUNT 

REAL,  I  NT  E  NT ( I  N) ,  Dl  ME  NS  I ON( VAL UE S COUNT )  ::  VALUES 

I  NTEGER,  I  NT  E  NT( I  NOUT)  :  :  ERRORNO 

CHARACTERf  256) ,  I  NT E NT ( I  NOUT )  ::  ERRORMSG 

!  V  a  r  i  a  b  I  e  s  to  hold  our  Para  me  ter  values 

REAL  WaterL,  Hb,  Cri  tWLev,  del  t  a  _  h ,  he,  Ga  mma  S ,  Ga  mma  W,  Di  mp,  MB  A 1  _  5  A I  I  I  R,  M  B  A 1  _  5  A I  I  I  S 

IVariables  to  hold  the  Indexes  of  our  parameters 
I  NTEGER 

I  X  _  Wa  t  e  r  L ,  I  X_Hb,  I  X_  C  r  i  tWLev,  I  X_  GammaS,  I  X_GammaW,  I  X_  D  i  mp,  I  X_  MB  A 1  _  5  A I  I  I  R,  I  X_  MBA1_  5AI  I  I  S 

IWARNING  -  Para  me  ter  names  below  may  change 
I  X  Wat  er  L=l  ndexOf  (  1  Wat  er  L1  ) 

I X'Hb=l  ndexOf  (  '  Hb'  ) 

I  X’Cr  i  t WL e v  =  I  ndexOf  (  1  Cr  I  t  WLev1  ) 

I  X_  Ga  mma  S  =1  n  d  ex  Of  (  '  Ga  mma  S '  ) 

I  X  _  G  a  mma  W=  I  n  d  e  x  Of  (  '  G  a  mma  W'  ) 

I  X_ D i  mp  =1  ndexOf  (  1  Di  mp1  ) 

I  X'MBAl  5  A I  I  I  R  =  l  ndexOf  (  '  MBA1  5AI  MR1) 

I  X”  MBA  1”  5  Al  I  I  S=l  ndexOf  (  '  MB  A 1  ~  5  Al  I  I  S'  ) 

I  Check  that  all  our  Para  me  ter  Names  have  been  recognised 
I  F  ( I  X  Wat  er  L<0.  OR.  I  X  Hb  <0 .  OR.  I  X  Cri  tWLev<0.  OR.  I  X  Ga  mma  S  <0 .  OR,  & 

'IX  GammaWcO.  OR'.  I  X  DI  mp  <0  ."OR .  IX  MB  Al  5  Al  I  I  R<0 .  OR.  I  X  MB  Al  5AI  I  I  S  <0 )  THEN 
E  RRORNO  =  l 

ERRORMSG=" One  or  more  Parameter  names  unrecognised" 

LSEBA1  5AM  I  =0.  0 
RETURN" 

END  I  F 

I  Get  the  parameter  values  that  we  need 
Wat  er  L  =VAL  UE  S  (  IX  Wat  er  L) 

Hb  =VAL UE S (  IX  Hb) ' 

Cri  t  WLev=VALUES(  I  X  Cri  tWLev) 

Ga  mma  S  =V A L UE  S ( I  X_  Ga  mma  S ) 

Ga  mma  W=VAL  UE  S  (  I  X  GammaW) 

Di  mp  =  VAL  UE  S (  I  X  DT  mp) 

MB Al  5  Al  I  I  R=VAL  UE  S ( I  X  MBA1  5AI  I  I  R) 

MB  AT  5  AM  I  S=VALUES(  I  X’  MB  A 1  ~  5AM  I  S) 

he  =(  Ga  mma  S  -  Ga  mma  W)  /  Ga  mma  W*  DI  mp 
del  t  a _ h =Wa  t  e r  L -  Hb 

LSEBAI5AI  I  I  =MBA1_  5AI  I  I  R  *  h  c - M  B  A 1 _  5  A I  I  I S  *  d  e I t  a  _  h 

END  FUNCTI  ON  LSEBA1_5AI  I  I _ 


Figure  2.3  An  example  FORTRAN  subroutine  for  a  soil  uplifting  LSE  is  provided  above 


It  is  anticipated  that  users  will  add  additional  LSEs.  Instructions  for  adding  additional  LSEs  are 
provided  in  Appendix  3. 


2.2  LSE  parameters  and  uncertainties 

Each  LSE  has  a  number  of  different  parameters  associated  with  it,  and  most  parameters  are  required 
within  more  than  one  LSE.  HR  Wallingford  (2008)  provides  details  a  list  and  description  of 
parameters  associated  with  the  FLOODsite  Task  4  report  (Allsop  et  al,  2007)  that  are  included  in  the 
pre-coded  LSEs  supplied  with  the  software.  HR  Wallingford  (2009)  also  provides  details  which 
parameters  are  used  by  which  LSE  (this  information  is  contained  within  the  FailureModePamm.csv 
file  which  is  available  from  the  RELIABLE  package). 
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The  user  is  required  to  specify  a  statistical  distribution  function  and  its  associated  parameter  values  for 
each  variable  within  each  LSE  used  in  the  construction  of  the  fault  tree.  Distributions  supported  by  the 
software  are: 

■  Normal 

■  Lognormal 

■  Gamma 

■  Rayleigh 

■  Weibull 

■  Bivariate  gamma 

It  is  also  possible  to  specify  constant  parameter  values  by  omitting  any  distribution  information  (i.e. 
only  a  single  value  is  entered). 

Whilst  the  user  is  required  to  specify  the  distributions,  example  information  on  distributions  is  stored 
within  the  excel  spreadsheet.  It  is  important  to  note  that  the  information  on  distribution  functions  and 
parameter  values  should  not  be  considered  “default  settings”  and  the  user  is  encouraged  to  gather 
evidence  for  the  structure  under  investigation,  to  support  the  selection  of  distribution  function  and 
parameters. 

2.3  Fault  Trees 

The  reliability  calculator  requires  the  user  the  opportunity  to  define  a  fault  tree  that  shows  the  logical 
links  between  the  LSEs  for  the  specific  structure  under  investigation.  Construction  of  fault  trees  is 
undertaken  in  freely  available  software  OpenFTA  (http://www.openfta.com/).  This  software  outputs  a 
text  file  (with  .fta  extension)  which  is  read  in  by  the  reliability  calculator.  Example  output  from 
OpenFTA  is  shown  below  (Figure  2.3).  An  example  of  the  fault  tree  construction  within  the 
OpenFTA  environment  is  shown  in  Figure  2.4. 


FailureMode.ped 

5  NULL  0 
0 

0 

M  NULL  1 

6  Breach 
O  NULL  3 
M  NULL  1 

33  Instability  due  to  anchor  failure 
H  NULL  2 
N  Cbl.2d  0 
M  NULL  1 


Figure  2.3  Text  file  output  from  OPENFTA 
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Figure  2.4  Example  fault  tree  for  sheet  pile  wall  analysis 

An  OpenFTA  database  file  (FailureMode.ped)  that  stores  the  list  of  (FLOODsite  Task  4)  failure  modes 
that  have  already  been  coded  is  provided  with  the  reliability  calculator,  to  facilitate  construction  of 
fault  trees  that  are  compatible  with  the  reliability  calculator. 

2.4  Interface 

The  reliability  calculator  interface  is  shown  in  Figure  2.5. 


5 


R  LIABLE 


Version  1.0 
OI-Jul-08 


IJob  n: 

TE2100 

User 

BPG 
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26-Jun-0a 

INPUTS 

Structure  Name  EBIS  -  _  . 

Calculate 

Reset 

Parameters 


Name 

CrestH 

DtssCoeff 

GrassRevQ 

Grav 

InsSlopeAng 

KdKf 

Mannrigs 

MBA11R 

MBA1_1S 

OverPercen 

StormD 


Distribution 

F  (Fixed) 

F  (Fixed) 

N  (Normal) 

F  (Fixed) 

F  (Fixed) 

N  (Normal) 

F  (Fixed) 

F  (Fixed) 

F  (Fixed) 

F  (Fixed) 

N  (Normal) 


Distribution  Parameters 
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i  W////zr/////// 

105  0  07V////////A 

981  V/////>/////// 
1749  ///////yy////// 

0  4  0  03  '////// S//S, 

1  V///////////// 
i  W/////r/////// 


8  'S/////J 
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HR  Wallingford  T  Delft 


Working  with  water 


ri.0  GSsite  u 


flood  risk 

management 


OUTPUTS 


Failure  Mode  Summary 
Name  Calls  Fal  Ratio 


1  Results 

Annual  ReliabiBty 

1 

Annual  Proba bitty  of  Falure  0 

Converged 

YES 

Convergence  Factor 

0 

Samples 

2200 

Time  Taken  (seconds) 

0 

Convergence 

control 

Min  Samples 

2000 

Max  Samples 

100000 

Interval 

200 

Convergence  Factor 

0.001 

Successive  Intervals 

1 

Control 

Events  Per  year 

1 

Optknise 

YES 

Log  each  Sample 

NO 

Log  each  LSE  call 

NO 

Failure  Mode  Summary 

Yes 

FRAGILITY  CURVE 

WaterL 

Probability 

This  interface  is  used  to  input  information  relating  to: 

■  The  name  of  the  structure  (determined  from  a  drop  down  list  which  uses  information 
from  OpenFTA  files  with  the  pre-constructed  fault  tree. 

■  Distribution  functions  for  each  parameter 

■  Parameters  for  the  distribution  functions 

■  The  number  of  samples  required  for  the  Monte-Carlo  simulation 

■  The  required  accuracy  of  the  calculation  (convergence) 

More  detailed  information  on  using  the  interface  and  relating  to  these  inputs  are  provided  in 
Appendix  5.  Information  on  extension  of  the  tool  to  include  additional  structure  fault  frees  and 
failure  modes  is  provided  in  Appendix  3. 

3.  Potential  improvements 

Potential  refinements  to  the  software  include  those  identified  in  Table  3. 1 


Table  3.1  Potential  refinements  to  RELIABLE 


Reliability  tool  component 

Improvement 

User  interface 

Provide  a  more  robust  user  interface 

Failure  mechanisms  and  dependencies 

■  Link  hydrodynamic  models  (through 

emulators) 

■  Link  to  Finite  element  geotechnical 
failure  mechanisms  (through  emulators) 

Fault  free 

■  Build  in  capability  rather  than  external 
software 

Statistical  models 

■  Include  correlation  and  extreme  values  in 
the  hydraulic  boundary  conditions 

■  Extend  the  range  of  distribution  functions 

■  Include  correlations  between  other 

random  variables 

■  Include  Autocorrelations  and  system 
effects 

Numerical  integration 

■  Extend  to  include  FORM/SORM  for 
example 

Sensitvity  analysis 

■  Variance  decomposion 
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ABSTRACT:  Embankments  are  constructed  for  the  retention  of  water  for  irrigation  and  supply,  and  for  pro¬ 
tecting  people,  land,  and  property  from  flooding.  Failure  of  any  embankment  poses  risks  to  people  and  prop¬ 
erty  nearby  and  the  services  provided  by  the  embankment.  The  review  of  breaching  of  embankments  in  this 
paper  identifies  significant  issues  in  the  parameterisation  of  the  processes  in  existing  models  and  the  data  used 
for  calibration.  This  paper  describes  the  development  of  a  new  model  the  failure  of  an  embankment  that  can 
simulate  breach  formation,  and  hence  consequent  risks,  more  reliably  than  existing  models.  The  model  uses 
the  standard  principles  of  hydraulics,  sediment  transport  and  soil  mechanics  and  introduces  a  new  methodol¬ 
ogy  to  model  the  lateral  growth  of  the  breach  based  upon  a  combination  of  continuous  erosion  and  mass  in¬ 
stability.  The  model  can  simulate  the  failure  of  different  embankments,  either  homogeneous  or  composite,  by 
overtopping  or  piping,  and  includes  a  probabilistic  distribution  for  simulating  embankment  condition  and  soil 
parameters.  The  model  has  been  tested  using  both  experimental  and  real  failure  data,  with  modelling  results 


showing  reasonable  agreement  with  observed  values  for 


1  INTRODUCTION 

Embankments  are  constructed  for  the  retention  of 
water  for  irrigation  and  supply,  and  the  protection  of 
people,  land,  and  property  from  flooding.  Failure  of 
any  embankment  poses  risks  to  people  and  property 
nearby  and  the  services  provided  by  the  embank¬ 
ment.  The  ability  to  maintain  assets,  and  provide  an 
acceptable  standard  of  service  for  water  supply  and 
flood  defence  therefore  depends  on  understanding 
and  predicting  performance  of  the  embankments  un¬ 
der  all  conditions.  Tools  currently  available  for 
simulating  embankment  failure  are  not  very  accurate 
(e.g.  Mohamed  et  al,  2001)  and  can  only  be  used  for 
indicative  assessments.  Consequently,  the  prediction 
of  flood  risk  from  embankment  breach  may  be 
similarly  inaccurate.  This  applies  particularly  in  the 
critical  zone  close  to  the  dam  or  embankment  where 
the  risk  to  life  is  greatest.  Flooding  from  the  failure 
of  the  Teton  dam  in  1976,  and  from  the  Mississippi 
and  Missouri  rivers,  in  1993,  and  more  recently,  the 
Yangtze  River  in  China  during  1998,  are  examples  of 
these  hazards. 

The  prediction  of  potential  breaches  and  the  con¬ 
sequent  flooding  are  thus  important  steps  in  manag¬ 
ing  the  risk  from  potential  embankment  failure. 
However,  breach  simulation  and  breach  parameter 
prediction  are  considered  to  contain  the  greatest  un¬ 


a  range  of  different  scenarios. 

certainty  of  all  aspects  of  dam  break  flood  forecast¬ 
ing  (Wurbs,  1987,  Singh,  1996,  and  Morris,  2000). 

A  research  project  to  investigate  breach  formation 
through  embankments  has  been  undertaken  at  HR 
Wallingford,  with  the  following  objectives: 

1  To  review  existing  methodologies  for  modelling 
embankment  breaching  processes. 

2  To  identify  any  weaknesses  within  these  method¬ 
ologies  and  determine  gaps  in  current  knowledge 
and  understanding. 

3  To  develop  a  new  methodology  that  improved  the 
accuracy  of  prediction  of  breach  formation 
through  embankment  dams  and  flood  embank¬ 
ments. 

This  paper  presents  an  overview  of  the  project 
giving  a  summary  of  weaknesses  identified  in  cur¬ 
rent  methodologies,  an  overview  of  a  new  method¬ 
ology  for  predicting  breach  formation,  and  an  as¬ 
sessment  of  model  performance  against  a  variety  of 
test  cases. 


2  STATE-OF-THE-ART  AND  FIMITATIONS 

In  the  last  forty  years,  many  models  have  been  de¬ 
veloped  to  simulate  the  embankment  breaching  pro¬ 
cess.  In  spite  of  this,  the  current  state-of-the-art  for 
predicting  the  breaching  process  still  has  many  un¬ 
certainties.  A  detailed  review  of  existing  models 


(Mohamed,  2002)  has  revealed  the  several  weak¬ 
nesses  and  gaps  within  the  modelling  process,  which 
are  discussed  in  the  sections  below. 

2. 1  Breach  Initiation 

Little  quantitative  information  is  known  about  the 
breach  initiation  processes  for  overtopping  or  piping 
failure.  Determining  how  the  breach  initiates  will 
help  in  reliably  determining  how  long  it  takes  for  a 
breach  to  develop  to  a  critical  point.  This  in  turn  can 
help  emergency  planners  in  establishing  flood  risk 
and  potential  warning  times  for  areas  downstream  of 
an  embankment. 

2.2  Breach  Location 

All  of  the  models  reviewed  (Mohamed,  2002)  as¬ 
sumed  a  breach  located  centrally  within  a  dam  or 
embankment.  However,  some  failure  cases  showed 
that  breaching  might  occur  near  an  abutment  rather 
than  in  the  middle  of  a  dam.  Examples  include  the 
failures  of  Teton  Dam  (Jansen,  1980),  Baldwin  Hills 
Dam  (Hamilton  et  al,  1971),  La-Josefina  (Abril, 
2001),  and  the  Euclides  da  Cunha  failure  (Hughes, 
1981).  The  breach  growth  and  hence  outflow  from  a 
centrally  located  breach  is  likely  to  be  different  from 
a  ‘side’  breach  in  terms  of  time  to  peak  discharge, 
peak  value,  and  hydrograph  shape. 

2.3  Data  For  Calibration  and  Verification 

Most  existing  models  were  calibrated  or  verified 
using  either  or  both  data  sets  from  the  Teton  Dam 
failure  or  the  Huaccoto  landslide  in  Peru.  The  docu¬ 
mented  data  for  these  two  events  is  not  very  de¬ 
tailed.  For  instance,  estimated  peak  outflow  from 
the  Teton  Dam  failure  ranged  from  45,000  to  80,000 
m3/s.  It  was  also  noticeable  that  some  authors  veri¬ 
fied  their  models  with  data  from  the  piping  failure  of 
the  Teton  Dam  in  spite  of  developing  their  models  to 
only  simulate  an  overtopping  failure.  Some  also 
calibrated  their  ‘central’  breach  models  (i.e.  unre¬ 
stricted  breach  growth)  with  side  breach  failure  data 
(i.e.  erosion  was  restricted  on  one  side  by  rock 
abutments).  These  wide  ranging  inconsistencies 
support  the  need  for  good  quality  data  sets  (such  as 
large-scale  experimental  data)  for  the  calibration  and 
verification  of  breach  models. 

2.4  Breach  Morphology 

Two  common  assumptions  in  many  of  the  existing 
models  are  a  constant  shape  (e.g.  rectangular,  trape¬ 
zoidal,  or  parabolic)  of  the  breach  and  the  uniform 
erosion  of  the  breach  section  during  the  formation 
process.  These  assumptions  simplify  the  equation(s) 
used  to  update  the  breach  section  at  each  time  step, 
but  they  seem  to  be  physically  unrealistic.  Assuming 


uniform  erosion  throughout  the  section  means  that 
the  part  of  the  breach  above  the  water  surface  will 
erode  at  the  same  rate  as  that  submerged,  which  is 
obviously  incorrect.  It  also  means  that  the  sides  of 
the  breach  below  the  water  surface  will  also  erode  at 
the  same  rate  as  the  breach  base,  which  is  inconsis¬ 
tent  with  the  flow  stress  distribution  along  the 
breach  sides.  Longitudinal  growth  of  the  breach  was 
assumed  to  be  parallel  to  the  downstream  face  in 
some  of  the  models  (BREACH  (Fread,  1988),  BRES 
(Vissser,  1998)  models).  This  representation  is  not 
compatible  with  the  assumptions  of  continuity  for 
sediment  transport  (Mohamed,  2002).  For  example, 
the  BREACH  model  computes  the  flow  depth  and 
velocity  along  the  downstream  face  of  the  embank¬ 
ment  using  the  steady  uniform  flow  equation.  This 
flow  condition,  if  combined  with  the  sediment  con¬ 
tinuity  equation,  will  not  give  parallel  retreat  of  the 
downstream  face. 

2.5  Hydraulics  of  the  Flow  over  the  Embankment 

Most  existing  breach  models  use  two  techniques  to 
simulate  flow  over  the  crest  and  on  the  downstream 
face  of  the  embankment.  These  are: 

-  the  broad  crested  weir  equation  and 

-  the  1-D  Saint  Venant  equations  (a  simplification) 
The  Saint  Venant  equations  incorporate  the  follow¬ 
ing  assumptions  (Cunge  et  al,  1980): 

-  The  flow  is  one-dimensional. 

-  The  water  pressure  is  hydrostatic. 

-  Boundary  friction  and  turbulence  effects  can  be 
accounted  by  steady  state  flow  resistance  laws. 

-  The  average  channel  bed  slope  is  small. 

It  is  clear  that  the  second  and  fourth  assumptions 
may  not  be  applicable  for  breaching  of  embank¬ 
ments.  Since  the  streamline  curvature  is  not  small, 
vertical  acceleration  may  not  be  negligible.  Also, 
the  downstream  face  of  the  embankment  may  be 
considered  as  steep  in  hydraulic  terms.  In  the  deri¬ 
vation  of  the  broad  crested  weir  formula,  the  curva¬ 
ture  of  the  flow  is  taken  into  consideration.  The 
weir  formula  is  thus  often  used  to  calculate  the  flow 
over  the  crest  since  it  accounts  for  the  acceleration 
of  the  flow  to  the  critical  point  on  the  crest. 

The  steady  non-uniform  flow  equations  have  also 
been  used  to  compute  the  water  depths,  velocities, 
and  energy  slope  on  the  downstream  slope  due  to  the 
short  reach  of  the  breach  channel  and  its  steep  slope 
and  their  relatively  simpler  computations  compared 
to  the  Saint  Venant  equations. 

2.6  Sediment  Transport  Equations 

The  selection  of  a  sediment  transport  equation  to  be 
used  in  any  mobile  bed  problem  is  difficult  and  is 
typically  based  on  professional  judgement,  previous 
experience,  or  even  personal  preference.  When  con- 


sidering  the  breaching  process,  the  problem  becomes 
even  more  difficult.  Most  existing  sediment  trans¬ 
port  equations  were  derived  for  steady  state  sub- 
critical  flow  conditions,  for  specific  types  of  sedi¬ 
ment,  and  for  a  certain  range  of  sediment  diameters. 
These  conditions  are  likely  to  be  violated  during  the 
breaching  process  since  conditions  are  typically  un¬ 
steady,  supercritical  flow,  and  with  a  wide  variety  of 
soil  types  used  for  embankment  construction.  Re¬ 
search  in  the  area  of  the  unsteady  non-uniform 
sediment  transport  is  still  in  its  early  stages  and  more 
work  is  required  in  order  to  achieve  reliable  results 
that  could  be  used  for  simulation  of  problems  such 
as  breach  formation.  However,  in  the  absence  of 
any  other  method  to  predict  the  sediment  transport, 
careful  selection  from  the  existing  sediment  trans¬ 
port  formulae  might  be  undertaken.  On  selection  of 
these  formulae,  the  following  might  be  taken  into 
consideration: 

-  Their  applicability  to  flow  on  steep  slopes  and  for 
supercritical  flow. 

-  Their  derivation  (e.g.  based  on  dam  breach  ex¬ 
perimental  data?). 

2.7  Geo-Mechanics  of  the  breach 

In  all  of  the  breaching  experiments  that  were  re¬ 
viewed  by  the  authors  during  this  research,  instabil¬ 
ity  of  the  breach  sides  was  observed  during  the 
breaching  process.  Most  existing  models  do  not 
consider  this  process  which  means  that  they  neglect 
a  process  that  is  likely  to  be  vital,  and  thus  the  cali¬ 
bration  and  validation  of  these  models  must  be 
questioned.  The  models  that  consider  this  process 
use  very  simplified  assumptions  e.g.  BREACH 
(Fread,  1988)  and  BEED  (Singh,  1997).  Assuming 
constant  breach  shape  and  uniform  erosion  of  this 
section  also  affects  the  accuracy  of  the  slope  stabil¬ 
ity  calculations  since  lateral  erosion  will  tend  to 
steepen  the  banks  and  the  breach  side  slope  will  get 
steeper  and  steeper  as  water  flows  through  the 
breach  (Osman  et  al,  1988).  This  means  that  the  side 
slope  of  the  breach  changes  throughout  the  simula¬ 
tion  and  hence  its  shape. 

2.8  Modelling  Composite  Embankments  and 
Surface  Protection  Layers: 

Despite  that  composite  embankments  represent  a 
significant  percentage  of  real  embankments  around 
the  world,  the  majority  of  existing  models  were  de¬ 
veloped  to  simulate  failure  of  homogeneous  em¬ 
bankments.  The  failure  of  composite  embankments 
might  involve  other  processes  such  as  core  wall  in¬ 
stability  and  mixed  sediment  transport  in  addition  to 
the  processes  encountered  with  homogeneous  em¬ 
bankments.  Moreover,  many  man  made  embank¬ 
ments  have  surface  protection  to  prevent  erosion  of 
the  embankment  faces  (Singh,  1996).  The  effect  of 


surface  protection  was  either  neglected  or  oversim¬ 
plified  in  many  of  the  existing  models  (BREACH 
(Fread,  1988)). 

2.9  Key  issues  for  the  research 

It  is  clear  that  there  are  many  gaps  in  our  knowledge 
for  reliably  predicting  breach  growth  and  location 
and  the  research  has  focused  on  the  following  issues: 

1  A  realistic  representation  of  breach  development 
during  the  breaching  process. 

2  A  more  accurate  analysis  of  the  breach  side  slope 
instability  process  and  the  transport  of  this  mate¬ 
rial  after  the  instability. 

3  A  methodology  to  model  the  failure  of  composite 
embankments,  including  the  effect  of  embank¬ 
ment  protection  layers  on  breach  development. 

3  DESCRIPTION  OF  A  NEW  METHODOLOGY 
FOR  BREACH  SIMULATION 

3.1  Overtopping  of  Homogeneous  Embankments 

Adjusting  the  breach  shape  is  a  crucial  process  in 
any  embankment-breach  model.  Several  methods 
have  been  used  in  existing  models  that  simulate 
breach  top  width  adjustments.  A  new  method  is  pro¬ 
posed  to  predict  top  width  adjustment.  The  process 
assumes  a  rectangular  initial  shape  of  the  breach,  as 
water  flows  into  the  breach  its  shape  and  side  slope 
will  change  as  shown  in  Figure  1(B).  The  bottom 
width  and  the  breach  depth  will  increase  as  the  water 
erodes  the  section  sides  and  bottom.  The  top  width 
will  not  change  significantly  and  can  be  assumed 
constant  until  slope  instability  is  encountered. 


Figure  1:  (A)  initial  breach  shape  (B)  Hypothetical  breach 
shape  after  three  successive  time  steps. 


The  process  is  a  combination  of  continuous  ero¬ 
sion  and  discrete  mass  failures  due  to  side  slope  in¬ 
stability.  Continuous  erosion  is  calculated  by  using 
a  sediment  transport  formula  to  quantify  the  volume 
of  the  sediment  transported.  Then  by  analysing  the 
effective  shear  stress  distribution1  for  the  breach 
section  the  new  breach  shape  can  be  obtained,  as 
erosion  can  be  assumed  to  be  proportional  to  the  ef¬ 
fective  shear  stress.  The  new  shape  of  the  breach 
may  be  approximated  as  shown  in  Figure  2. 


Figure  2:  Approximated  shape  of  the  breach 


The  breach  section  is  updated  at  each  time  step 
assuming  that  maximum  lateral  erosion  (Db)  is  equal 
to  the  vertical  erosion  (Dh)  and  it  is  assumed  to  oc¬ 
cur  very  near  to  the  bottom  level  of  the  breach.  The 
top  width  is  kept  constant  until  slope  instability  is 
encountered.  The  stability  of  the  breach  side  slope  is 
analysed  by  taking  into  consideration  the  forces 
acting  on  the  slope,  and  variation  of  the  soil  density. 
A  factor  of  stability  (FOS)  is  obtained  using  the  fol¬ 
lowing  equation: 

FOS  =  (Stabilising  Forces  /  Destabilising  Forces)  (1) 
Where: 

The  stabilising  forces  are: 

-  Water  pressure  forces  in  the  breach  channel. 

-  Friction  forces. 

-  Cohesion  forces  (if  any). 

The  destabilising  forces  are: 

-  Gravity  forces. 

-  Pore  water  pressure  forces  in  the  embankment. 

The  nearly  vertical  sides  of  the  breach  (as  ob¬ 
served  in  both  lab  experiments  and  real  failures) 
suggest  that  slope  instability  failure  modes  might  be 
either  through  shear  or  bending  failure  (Mohamed, 
2002).  Both  these  modes  of  failure  can  lead  to  a 


near  vertical  failure  plane.  In  the  following  two  sec¬ 
tions,  a  description  of  each  failure  mode  is  given. 

3.1.1  Bending  Failure 

An  initial  rectangular  notch  is  assumed  on  the  crest 
and  the  downstream  face  of  the  embankment.  Water 
flows  through  the  initial  notch.  Flowing  water 
erodes  the  breach  sides  below  the  water  surface  and 
the  bottom  of  the  notch  and  undermines  the  slope. 
The  erosion  process  continues  until  slope  instability 
is  encountered.  A  tension  crack  develops  progres¬ 
sively  as  the  actual  tension  stress  exceeds  the  soil 
tension  strength.  The  soil  block  rotates  and  falls  into 
the  flowing  water.  Water  erodes  the  slumped  mate¬ 
rial  and  the  process  continues  until  the  reservoir  is 
depleted  or  the  breach  reaches  its  maximum  dimen¬ 
sions.  This  mode  of  failure  is  likely  to  occur  in  cohe¬ 
sive  embankments. 

The  following  assumptions  have  been  made  in 
developing  the  analysis  below  (Figure  3): 

-  Suction  is  neglected  in  the  zone  above  the  water 
level  within  the  embankment.  This  zone  is  con¬ 
sidered  dry. 

-  Changes  in  water  level  inside  the  embankment 
during  the  embankment  failure  time  are  small  and 
can  be  neglected. 


Moment(Mo)  =  W.e+ Ws.es  +Wu.eu  +H2.e2  -  Hj.e,  (2) 


Where: 

W 

:  Weight  of  dry  block  of  the  soil. 

Ws 

:  Weight  of  saturated  block  of  the  soil. 

Wu 

:  Weight  of  submerged  block  of  the  soil. 

H, 

:  Hydrostatic  pressure  force  in  the  breach 
channel 

h2 

:  Hydrostatic  pressure  force  inside  the  em¬ 
bankment. 

c,  es,eu 

:  Weight  forces  eccentricities. 

ei,  e2 

:  Hydrostatic  pressure  forces  eccentricities. 

L 

:  Length  of  the  failure  plane. 

1  Effective  shear  stress  equals  the  difference  between  the  total  shear 
stress  and  the  critical  shear  stress. 


Based  on  the  above  analysis,  the  maximum  actual 
tensile  stress  (at(actuai) )  on  the  plane  of  failure  can  be 
computed  as  follows: 

crl(actual)=(H2-H,)/L  +  6M0/L2  (3) 

Assuming  that  the  allowable  soil  tensile  strength 
(Ot(soii))  is  known,  then  ot(actuai)  is  compared  with 
at(soii)  and  if  (at(actuai)  >  at(Soii) )  then  failure  occurs. 

3.1.2  Shear  Failure 

A  similar  process  to  that  explained  above  occurs  for 
this  failure  mode,  however  the  slope  fails  due  to  ex¬ 
ceeding  the  shear  strength  of  the  soil  (Figure  4).  This 
mode  of  failure  is  likely  to  occur  in  a  non-cohesive 
embankment. 


Figure  4:  Forces  for  shear  failure. 


Making  similar  assumptions  to  the  bending  fail¬ 
ure  analysis  above: 


FOS 


c  *  L  +  H ,  tan  (f> 

W  +  Ws  +  Wu  +H2  tan 


(4) 


Where: 

c  :  Soil  cohesion. 

4>  :  Soil  angle  of  friction. 


3.1.3  Dealing  with  uncertainty  in  soil  properties 
and  construction  quality 

A  probabilistic  approach  is  used  to  take  into  account 
uncertainties  in  soil  properties  and  the  quality  of 
construction.  A  Sigmoid  function  has  been  used  to 
represent  a  probability  distribution  for  the  factor  of 
safety.  This  was  selected  since  it  allows  both  ex¬ 
tremes  of  100%  and  0%  to  be  represented  along  with 
various  ranges  of  distribution  in  between: 


f  (m  ) 


1 

1  +  e  a(m  _1) 


(5) 


The  value  of  m  represents  the  factor  of  safety. 
The  uncertainty  coefficient,  a,  controls  the  probabil¬ 
ity  distribution  and  may  represent  the  quality  or 
knowledge  of  materials  within  and  construction  of 
the  embankment  (e.g.  very  good,  good  or  poor  mate¬ 
rial  and  construction).  Three  different  uncertainty 


distributions  (Figure  5)  were  used  to  represent  the 
quality  of  materials  and  construction,  however,  other 
distributions  may  also  be  used  to  represent  varying 
conditions  and  uncertainty. 


Figure  5:  Probability  distribution  functions  for  the  slope 
stability  factor  of  safety 

3.1.4  Cohesive  Embankments: 

The  methodology  discussed  so  far  is  mainly  for  non- 
cohesive  embankments  with  some  apparent  or  con¬ 
ventional  cohesion.  For  cohesive  embankments,  the 
failure  process  might  be  different.  Hughes  (1981), 
Al-Qaser  (1991),  and  Hanson  (2000)  conducted 
laboratory  and  field  experiments  to  determine  the 
failure  process  of  cohesive  embankments  due  to 
overtopping.  They  observed  the  formation  of  an  over 
fall  or  steps  that  progressively  advanced  towards  the 
upstream  face  (Figure  6).  Hanson  (2000)  concluded 
that  the  erosion  process  and  soil  type  have  a  signifi¬ 
cant  effect  on  the  timing  and  rate  of  discharge  during 
overtopping  events.  The  observed  processes  were 
described  as  follows: 

-  Initial  downstream  surface  erosion. 

-  This  initial  erosion  progresses  into  stair-stepped 
multiple  over  falls. 

-  Over  falls  then  merge  into  a  single  upstream  mi¬ 
grating  head  cut. 

-  The  head  cut  then  migrates  upstream,  lowering 
the  crest  by  advancing  into  the  upstream  em¬ 
bankment  face. 


Figure  6:  Headcut  advance  mechanism  (Flanson  et  at,  2000) 
Of  all  the  models  reviewed  within  this  research, 
only  the  SITES  model  (Wahl,  1998)  simulates  the 


first  three  processes  described  above.  It  does  not, 
however,  model  the  fourth  process,  which  is  critical 
for  predicting  the  flow  of  water  from  a  breach. 

3.2  Overtopping  of  Composite  Embankments 

The  failure  of  composite  embankments  differs  from 
that  of  homogeneous  embankments,  because  of  the 
existence  of  less  erosive  layers  (such  as  a  clay  core) 
within  the  dam  body.  The  erosion  of  the  material 
behind  the  core  may  affect  the  stability  of  the  core 
and  could  eventually  lead  to  its  failure.  The  likely 
failure  mechanisms  of  the  core  wall  of  a  composite 
embankment  dam,  assuming  that  a  large  part  of  the 
downstream  face  has  been  eroded,  are: 

-  Sliding  of  the  clay  core  wall. 

-  Overturning  of  the  clay  core  wall. 

-  Bending  of  the  core  wall. 

In  the  following  sections  a  description  of  the  pos¬ 
sible  modes  of  failure  of  the  core  are  presented. 

3.2.1  Sliding  Failure: 

Figure  7  shows  the  forces  acting  on  the  core  after  a 
large  amount  of  the  downstream  body  material  has 
been  eroded.  These  forces  include: 

-  Active  soil  earth  pressure  forces  from  the  material 
behind  the  core  (FI). 

-  Water  pressure  forces  (F2). 

-  Weight  of  the  clay  core  material  above  the  failure 
plane  (F3). 

-  Weight  of  the  material  above  the  upstream  face  of 
the  clay  core  and  the  failure  plane  (F4). 

-  Weight  of  the  water  above  the  core  (F5). 

-  Earth  pressure  forces  on  the  core  sides  (F6). 


Figure  7:  Forces  acting  on  the  clay  core  wall. 

The  failure  plane  is  assumed  to  be  (as  shown  in 
Figure  7)  just  above  the  non-eroded  material  of  the 
downstream  face.  Forces  (FI)  and  (F2)  are  the  de¬ 
stabilising  forces,  while  the  forces  (F3),  (F4),  (F5), 
and  (F6)  are  stabilising  forces.  The  latter  forces  mo¬ 
bilise  the  friction  on  the  bottom  and  the  sides  of  the 
failure  plane.  The  cohesion  on  the  right,  left,  and 
bottom  sides  of  the  failure  plane  will  also  resist  the 
destabilising  forces.  The  ratio  of  the  stabilising  to 
the  destabilising  forces  is  calculated  and  if  it  is  less 
than  unity  the  core  is  considered  to  have  failed. 


3.2.2  Overturning  Failure: 

A  similar  analysis  to  that  for  sliding  failure  is 
used  here.  However,  the  stabilising  and  destabilising 
moments  are  compared.  If  the  ratio  between  them  is 
less  than  unity  the  core  is  considered  to  fail  by 
overturning. 

3.2.3  Bending  Failure: 

Hughes  (1981)  showed  that  the  core  failure  mecha¬ 
nism  could  be  as  shown  in  Figure  8.  He  supported 
this  by  the  observations  of  Jayasinghe  (1978)  who 
noticed  similar  cracks  to  those  shown  in  Figure  8  be¬ 
fore  failure  of  the  core  wall. 


restrained 

sides 


Figure  8:  Typical  cracking  of  a  clay  core  restrained  on  three 
edges  and  subject  to  overtopping  (Flughes  (1981)). 

Analysis  of  this  failure  mechanism  is  difficult  be¬ 
cause: 

1.  The  geometry  of  the  core  and  the  variation  of 
forces  with  the  height  of  the  core  wall. 

2.  The  difficulty  in  defining  the  support  type  on  the 
sides  and  the  bottom  of  the  core  wall  (hinged, 
fixed,  or  something  in  between). 

3.  The  difficulty  in  computing  the  load  distribution 
factors  of  the  core  sides  and  bottom. 

The  following  simplified  analysis  is  proposed  for 
this  problem.  The  clay  core  might  fail  by  a  bending 
failure  as  shown  in  Figure  9  if  the  left,  right,  bottom 
sides  can  be  considered  as  hinges.  Figure  9  (A) 
shows  the  failure  plane  in  that  case.  Figure  9  (B) 
shows  an  assumed  failure  plane  neglecting  the  hinge 
at  the  bottom  side.  This  assumption  can  be  justified 
since  the  breach  depth  is  likely  to  be  greater  than  the 
breach  width  at  that  stage  of  the  embankment  failure 
so  the  effect  of  the  hinge  at  bottom  can  be  neglected. 


Figure  10  shows  the  forces  acting  on  the  clay  core 
in  the  case  of  bending  failure  after  a  large  part  of  the 
downstream  face  material  has  been  eroded,  namely: 


-  Active  earth  pressure  from  the  material  behind 
the  core  (FI). 

-  Water  pressure  forces  (F2). 


Figure  10:  Forces  acting  on  the  core  wall  in  bending  failure. 

The  core  is  assumed  to  be  as  a  simple  beam  and 
the  max  bending  moment  and  stress  is  computed. 
The  bending  stress  can  be  compared  to  the  tension 
strength  of  the  clay  core  material.  If  the  bending 
stress  exceeds  the  tension  strength  of  the  core  mate¬ 
rial  then  failure  will  take  place. 

3.3  Piping  of  Embankments 

Piping  starts  when  water  flows  through  an  embank¬ 
ment  body  via  cracks,  cavities,  weak  layers  etc.  As 
flow  increases,  the  water  starts  to  remove  soil  parti¬ 
cles  until  a  pipe  (tunnel)  is  formed  through  the  body 
of  the  embankment  ( 

Figure  11  (A)). 


After  the  formation  of  a  pipe  inside  an  embank¬ 
ment,  the  failure  process  can  be  divided  into  the 

following  processes: 

-  Erosion  of  material  in  the  pipe. 

-  Slumping  of  the  downstream  face  above  the  pipe 
(Figure  11  (B  and  C)). 

-  Collapse  of  the  embankment  crest,  under  its  own 
weight  or  by  water  pressure  Figure  1 1  ((D)). 

-  Erosion  of  the  embankment  body  in  a  similar 
manner  to  that  for  overtopping  breach  formation 
(Sections  3.1  and  3.2). 

-  In  the  following  sections  a  description  of  the 
modelling  procedure  for  each  process  is  given. 


3.3.1  Erosion  of  the  material  in  the  pipe. 

Figure  12  shows  a  typical  hydrograph  for  the  piping 
failure  of  an  embankment.  It  can  be  noted  that  dif¬ 
ferent  time  scales  are  associated  with  each  process. 
For  the  Teton  Dam  these  times  were  approximately 
as  follows: 

-  The  initiation  stage  was  at  least  2  days. 

-  Erosion  of  the  pipe  stage  until  the  collapse  of  the 
top  of  the  dam  was  about  4  hours. 

-  Emptying  of  the  reservoir  was  about  2  hours. 


The  initiation  stage  was  not  studied  in  this  research 
but  attention  was  focussed  on  growth  of  an  initial 
pipe. 


Figure  12:  Time  scales  for  piping  processes 

An  initial  pipe  is  assumed  to  establish  within  the 
dam  body  at  the  start  of  the  simulation.  The  flow 
through  the  pipe  is  calculated  using  the  orifice  flow 
equation  as  follows: 


= 

A  -y/2  g  (H  -  H  p)/hL 

(6) 

Where: 

Qb 

:  Flow  through  the  pipe. 

g 

:  Acceleration  due  to  gravity. 

A 

:  Pipe  cross  section  area. 

H 

:  Water  level  in  the  dam. 

HP 

:  Pipe  centre  line  elevation. 

hL 

:  Losses  due  to  friction  and  contraction. 

K  = 

o 

o 

+ 

b  'pj 

(7) 

Where: 

f  :  Friction  coefficient  determined  as  function  of  D50. 
L  :  Pipe  Length. 

D  :  Pipe  Diameter. 

The  factor  (0.05)  in  the  above  equation  represents 
the  coefficient  of  the  contraction  losses.  The  en¬ 
largement  of  the  pipe  was  modelled  using  the  sim¬ 
plified  procedure  used  by  Fread  (1988)  and  Paquier 
et  al  (1998).  The  pipe  is  enlarged  uniformly  along 
its  length  based  on  the  volume  of  sediment  eroded 


(Vs)  within  a  specified  time  step.  Vs  can  be  com¬ 
puted  as  shown  below: 

Vs  =  Q  sAt  (8) 

Where: 

Qs  :  Sediment  transport  rate. 

At  :  Time  step. 

It  should  be  noted  that  this  procedure  has  been 
used  in  the  absence  of  any  more  detailed  approaches. 
This  highlights  the  need  for  research  on  the  initiation 
of  piping.  This  simplified  procedure  can  be  used  for 
homogeneous  and  composite  embankment  dams. 

3.3.2  Slumping  of  downstream  face  material  above 
the  pipe 

As  the  pipe  develops,  material  from  the  downstream 
face  falls  into  the  pipe  and  is  swept  away  in  the 
flow.  This  mechanism  was  observed  during  the 
piping  failure  of  the  Teton  dam  in  1976.  Justin  et  al 
(1945)  and  Hagerty  (1991)  also  reported  this  proc¬ 
ess.  It  is  anticipated  that  the  material  falls  under  its 
own  weight.  The  vertical  failure  planes  observed 
during  the  Teton  Dam  failure  and  reported  by 
Hagerty  (1991)  suggest  that  this  is  likely  to  be  a 
shear  failure.  The  ratio  between  the  weight  of  the 
material  of  the  wedge  and  the  soil  strength  is  com¬ 
puted.  If  its  value  is  below  unity  then  the  wedge  of 
soil  above  the  top  of  the  pipe  is  assumed  to  be  un¬ 
stable. 

3.3.3  The  collapse  of  the  top  part  of  the  dam,  above 
the  pipe. 

As  material  from  the  downstream  face  slumps  into 
the  pipe,  the  width  of  material  above  the  pipe  re¬ 
duces.  If  water  forces  are  high  enough  to  exceed  the 
shear  strength  of  the  embankment  material  above  the 
pipe,  then  this  ‘wedge’  of  material  may  fail.  The 
‘wedge’  may  also  collapse  under  it  is  own  weight. 
This  mechanism  was  observed  during  failure  of  the 
Teton  Dam.  For  more  information  about  these  fail¬ 
ure  mechanisms,  the  reader  is  referred  to  Mohamed 
(2002). 

3.3.4  Erosion  of  the  dam  body. 

The  process  here  is  similar  to  that  described  in  Sec¬ 
tions  3.1  and  3.2  of  this  paper,  which  is  a  combina¬ 
tion  of  erosion  and  slumping  of  material. 


4  MODEL  PERFORMANCE: 

Three  sets  of  data  were  selected  to  test  different  as¬ 
pects  of  the  proposed  new  modelling  methodology. 
A  detailed  description  of  each  case  is  given  in  Mo¬ 
hamed  (2002)  these  are  briefly  summarised  below. 
The  test  data  was  obtained  from: 

-  CAD  AM  project  Munich  proceedings  (1998). 


-  The  Teton  Dam  failure  (Fread  (1988),  and  Jansen 

(1980)). 

4.1  Case  1 

To  check  the  model  performance  for  modelling  the 
overtopping  failure  of  homogeneous  embankments, 
a  set  of  data  was  used  from  a  series  of  laboratory  ex¬ 
periments  that  were  carried  out  at  the  Federal  Armed 
Forces  University  in  Munich  (Bechteler  et  al,  1998) 
to  explore  the  3D-development  of  a  dam  breach. 

Figure  13  shows  the  results  of  test  case  1.  Three 
simulations  have  been  carried  out  using  different 
sediment  transport  equations  (Visser  (1998),  Yang 
(1979),  and  Chen  and  Anderson  (1986)).  The  avail¬ 
able  measured  outflow  hydrograph  was  compared 
against  the  corresponding  simulation  results.  For 
Visser’ s  and  Yang’s  equations,  it  can  be  seen  that 
the  model  results  are  reasonably  in  agreement  with 
the  measured  results  in  terms  of  peak  outflow  value 
and  timing.  Chen’s  sediment  transport  equation  gave 
poorer  results  than  the  other  two.  The  jagged  nature 
of  the  plots  arises  from  breach  growth  through  bank 
instability. 
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Figure  13:  Test  case  1  outflow  comparison. 

Within  the  CAD  AM  project,  several  models  were 
tested  against  the  data.  A  comparison  between  these 
models  and  the  developed  model  is  shown  in  Figure 
14  below. 


Apart  from  the  results  of  Broich’s  model  (which 


was  calibrated  against  the  test  case  data)  the  model 
performance  is  generally  better  than  the  other  mod¬ 
els  in  terms  of  predicted  peak  outflow  time  and 
value. 

4.2  Case  2 

To  check  the  model  performance  in  modelling  the 
overtopping  failure  of  composite  embankments,  the 
data  was  used  from  a  fuse  plug  breach  test  that  was 
conducted  in  1982  in  the  spillway  chute  of  the  Ya- 
hekou  Dam  in  China.  The  fuse  plug  embankment 
was  made  of  sand  fill  with  a  clay  core,  5.6  m  high 
and  41  m  long  at  the  top  and  31m  long  at  the  bot¬ 
tom,  with  a  crest  width  of  4m,  creating  a  reservoir 
with  a  capacity  of  46,000  m3. 

Two  simulations  were  undertaken  using  this  test 
case  data  to  investigate  the  effect  of  core  failure 
mechanisms.  The  first  considered  only  sliding  and 
overturning  failure  modes  (A)  whilst  the  second  in¬ 
cluded  the  bending  failure  mode  (B)  as  well.  This 
was  undertaken  to  see  whether  the  bending  failure 
mechanism  is  a  critical  mode  of  failure  or  not,  since 
modelling  this  failure  was  based  on  a  simplification 
of  the  real  processes  (See  Figure  15). 


Figure  15:  Test  case  4  outflow  comparison. 


The  predicted  peak  flow  and  its  time  for  simula¬ 
tion  A  was  closer  to  the  measured  value  than  those 
predicted  in  simulation  B.  However,  the  predicted 
time  of  the  first  failure  of  the  core  was  earlier  than 
what  was  observed  during  the  experiment  for  the 
two  simulations.  For  simulation  A,  the  magnitude  of 
the  flow  at  the  first  core  failure  was  significantly 
higher  than  the  measured  value.  Whilst,  for  simula¬ 
tion  B,  it  was  slightly  higher  than  the  measured 
value.  The  results  obtained  here  showed  that  further 
investigations  and  tests  are  needed  to  check  the  per¬ 
formance  of  the  model  of  the  composite  embank¬ 
ments. 

4.3  Case  3 

The  data  of  Teton  dam  failure  was  used  to  test  the 
model  performance  in  modelling  piping  failure  of 
embankments.  The  Teton  dam  failed  on  the  6th  of 


June  1976.  The  dam  was  a  91.45  m  high  earth  dam 
with  a  914.5  m  crest  length,  10.65  m  crest  width, 
and  79.85  m  depth  storing  about  307  million  m3  of 
water  at  the  time  of  failure  (Fread,  1988).  Failure 
occurred  close  to  the  right  rock  abutment,  which 
potentially  restricted  lateral  growth  of  the  breach. 

Modelling  the  piping  failure  was  undertaken  as¬ 
suming  a  small  initial  pipe  within  the  dam,  and  then 
a  central  breach  or  a  side  breach  after  the  top  of  the 
dam  collapses.  Figure  16  shows  the  results  obtained 
from  both  simulations  using  Yang’s  equation.  It  can 
be  seen  that  the  central  breach  simulation  gave  a 
higher  peak  outflow  value  than  the  side  breach  case. 
The  results  of  the  central  and  side  breach  simula¬ 
tions  highlight  the  importance  that  the  breach  loca¬ 
tion  and  the  valley  shape  might  have  on  the  peak 
outflow  value.  This  also  raises  questions  as  to  the 
validity  of  models  (such  as  BREACH,  BEED)  which 
were  calibrated  using  the  Teton  dam  failure  data  yet 
only  assumed  a  centrally  located  breach. 


Figure  16:  Test  case  3  outflow  comparison. 


To  test  the  effect  of  using  different  sediment 
transport  equations  on  the  model  results,  the  Visser, 
Yang,  and  Chen  and  Anderson  equations  were  also 
used  assuming  a  side  breach  failure.  Figure  17 
shows  the  results  obtained  using  these  three  equa¬ 
tions. 
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Figure  17:  Test  case  5  outflow  comparison  using  different 
sediment  equations. 

The  peak  outflow  values  obtained  using  Visser 
and  Yang  equations  were  within  the  estimated  range 


of  the  peak  outflow  value.  However  the  predicted 
value  using  Yang’s  equation  was  less  than  that  ob¬ 
tained  using  Visser’s  equation.  The  value  obtained 
using  Chen  and  Anderson’s  equation  was  below  the 
estimated  range.  It  can  also  be  seen  that  the  time  to 
peak  obtained  using  each  equation  varied  signifi¬ 
cantly.  This  variation  highlights  the  importance  of 
selecting  a  sediment  transport  equation  in  dam 
breach  problem  and  their  effect  on  the  results. 


5  CONCLUSIONS: 

There  are  several  weaknesses  in  the  commonly-used 
methods  for  simulating  breaching  of  embankments. 
These  arise  from  what  now  appear  to  be  unrealistic 
simplifications  of  the  breach  growth  processes  and 
from  using  as  calibration  data  information  from  past 
failures  where  the  valley  topography  or  failure  mode 
were  different  from  those  included  in  the  model. 
Thus  the  situation  is  that  there  remains  a  high  degree 
of  uncertainty  on  the  assessed  outflow  hydrograph  of 
the  water  impounded  behind  an  embankment  during 
a  failure.  This  has  important  social  and  economic 
impacts,  particularly  for  risk  management  and  emer¬ 
gency  planning. 

A  new  methodology  has  been  put  forward  in  this 
paper  based  upon  well-established  principles  of  soil 
and  fluid  mechanics,  although  uncertainty  still  re¬ 
mains  in  some  aspects  of  the  parameterisation  of  the 
model.  The  evidence  from  the  documentation  of 
actual  failures  is  that  the  failure  of  an  embankment 
progresses  through  general  erosion  of  the  embank¬ 
ment  material  and  mass  failures  caused  by  soil  slope 
instability.  The  novelties  in  the  model  described  in 
this  paper  includes  updating  the  shape  of  the  eroding 
part  of  the  breach  at  each  time  step  rather  than  as¬ 
suming  a  fixed  shape  and  allowing  for  uncertainty  in 
the  knowledge  of  soil  properties  through  a  probabil¬ 
istic  approach.  The  model  has  not  been  calibrated 
on  a  particular  data  set  but  has  been  shown  to  per¬ 
form  well  on  a  range  of  data  obtained  from  labora¬ 
tory  investigation,  field  scale  measurements  and 
historic  failures. 

Whilst  the  model  presented  here  represents  and 
advance  on  other  available  methods,  several  issues 
remain  for  further  research  including: 

-  the  failure  processes  of  embankments  with  sub¬ 
stantial  flow  parallel  to  the  crest 

-  the  initiation  of  piping  failures 

-  fracturing  mechanisms  of  exposed  clay  cores 

-  the  influence  of  downstream  topography  on  the 
failure  process 

-  the  sediment  transport  rates  for  highly  unsteady, 
non-uniform  flows,  and 

-  sediment  transport  for  natural  mixtures  arising 
from  landslide  dams. 
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This  paper  describes  a  rapid  embankment  breach  modelling  method  that  has  been  developed  as  part  of  the 
FRMRC2  programme  to  improve  current  methods  for  predicting  breach  outflow  volumes  and  hence  to  improve 
the  accuracy  of  flood  risk  assessments.  Application  of  the  method  led  to  the  development  of  a  fast  model 
called  AREBA  which  rapidly  predicts  the  breach  hydrograph  of  homogeneous  embankment  failures  whereby 
discharge  through  the  breach  is  controlled  by  the  breach  dimensions.  AREBA’ s  surface  erosion  failure  mode, 
and  piping  failure  mode  have  been  validated  against  the  IMPACT  project  data  (www.impact-project.net),  and 
benchmarked  against  HR  BREACH  model  version  4.1.  Validation  showed  that  AREBA  is  able  to  predict 
the  breach  hydrograph  within  the  bounds  of  uncertainty  that  originate  from  uncertainty  in  the  model  input 
parameters. 


1  INTRODUCTION 

System  risk  models  for  assessing  flood  risk  at  a  catch¬ 
ment  or  larger  scale  are  used  to  determine  how  in¬ 
vestment  in  flood  control  measures  might  be  opti¬ 
mised.  The  system  risk  approach  requires  a  rapid  and 
accurate  method  of  assessing  the  impact  of  a  possi¬ 
ble  flood  event  on  a  given  system  of  flood  protec¬ 
tion.  One  part  of  that  process  is  the  prediction  of  a 
flood  volume  through  an  embankment  breach.  The 
AREBA  model,  developed  as  part  of  the  FRMRC2 
programme,  rapidly  predicts  a  breach  hydrograph,  fi¬ 
nal  breach  width  and  final  breach  depth  for  grass  cov¬ 
ered  embankments  that  fail  due  to  overflow  or  pip¬ 
ing.  Within  the  overflow  failures,  surface  erosion  fail¬ 
ures  and  headcut  erosion  failures  are  distinguished. 
The  AREBA  model  is  limited  to  predicting  failures 
of  homogeneous,  trapezoidal  shaped  flood  embank¬ 
ments  and  requires  soil  parameters,  an  embankment 
geometry,  and  upstream  and  downstream  boundary 
conditions,  as  input.  A  box  shaped  reservoir  with  an 
inflow  and  outflow  is  assumed  to  allow  the  use  of  a 
triangular  or  trapezoidal  hydrograph  as  the  upstream 
boundary  condition  for  the  breach  model.  The  down¬ 


stream  boundary  condition  consists  of  a  user-defined 
flood  area  which  for  simplicity  has  been  assumed  con¬ 
stant  over  the  height  of  the  embankment.  AREBA  can 
however  be  easily  adapted  to  deal  with  more  compli¬ 
cated  boundary  conditions  like  stage-area  curves  or 
specific  river  hydrographs.  This  paper  addresses  the 
constitutive  equations  and  validation  of  the  surface 
erosion  and  piping  failure  modes  and  shortly  consid¬ 
ers  the  headcut  failure  mode.  The  paper  also  discusses 
the  accuracy  of  AREBA  with  respect  to  the  assump¬ 
tions  made.  Section  2  focuses  on  the  simplifications 
made  with  respect  to  the  breach  physics  to  obtain  a 
fast  model.  The  test  results  of  AREBA  are  presented 
in  Section  3  and  discussed  in  Section  4. 

2  MODELLING  METHODOLOGY 

2.1  Overflow  failures 

In  the  case  of  overflow  failure  the  overall  erosion  be¬ 
haviour  is  a  balance  between  the  rate  at  which  the 
crest  erodes  downward  due  to  flow  over  the  crest,  and 
the  rate  with  which  the  landside  slope  erodes  towards 
the  waterside  slope  due  to  the  flow  over  the  landside 
slope.  Within  overflow  failures  AREBA  distinguishes 
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between  headcut  erosion  and  surface  erosion  (see  fig¬ 
ure  1).  In  the  case  of  headcut  erosion,  the  downward 


Surface  erosion  Headcut  erosion 

Figure  1:  Illustration  of  surface  and  headcut  erosion 

erosion  rate  of  the  crest  is  negligible  with  respect  to 
the  rate  at  which  the  landside  slope  retreats  towards 
the  waterside  slope  (Temple  et  al.,  2005).  In  the  case 
of  surface  erosion,  both  processes  play  an  equally  im¬ 
portant  role.  In  AREBA,  the  choice  of  failure  mode 
is  user  defined.  A  broad  assumption  would  be  to  as¬ 
sume  headcut  failure  for  a  strong  erosion  resistant  soil 
such  as  a  strong  clay,  and  surface  erosion  for  a  weaker, 
more  erodible  soil  with  higher  sand  or  gravel  con¬ 
tent.  However  the  point  of  transition  between  these 
two  processes  as  a  function  of  soil  erodibility  is  un¬ 
clear,  hence  where  uncertain,  both  failure  modes  may 
be  considered.  As  a  starting  point  to  calculating  the 
erosion  rates,  AREBA  calculates  the  discharge  over 
the  crest  using  the  following  broad  crested  weir  equa¬ 
tion  to  account  for  the  effects  of  non-hydrostatic  pres¬ 
sures  over  the  crest  due  to  flow  contraction. 

Qb  =  cwbh2\/2g(H  -  h2),  (1) 

where  Qb  is  the  discharge  over  the  crest  in  m3/s, 
cw  the  weir  coefficient,  H  —  h2  (m)  is  the  difference 
in  water  level  in  the  reservoir  H  and  the  waterlevel 
over  the  hydraulic  control  h2,  and  b  is  the  flow  width 
which  is  assumed  to  be  spatially  constant  (Nortier  and 
De  Koning,  1991).  The  hydraulic  control  is  assumed 
to  be  located  where  the  crest  meets  the  waterside 
slope.  The  water  depth  d  over  the  crest  is  expected 
to  be  equal  to  2/3 H.  h2  is  assumed  to  be  equal  to 
the  depth  d,  for  non-drowned  weir  flow,  and  equal  to 
the  downstream  water  level  for  a  drowned  weir  flow. 
Every  time  step,  the  upstream  and  downstream  water 
levels  are  updated  assuming  an  equal  drop  or  rise  in 
water  level  everywhere  in  the  reservoir  or  flood  plain. 
The  weir  coefficient  cw  is  a  combination  of  three  coef¬ 
ficients.  Before  the  landside  slope  has  reached  the  wa¬ 
terside  slope  cw  =  cw i  whereby  cwl  accounts  for  the 
effects  of  vertical  flow  contraction.  When  the  landside 
slope  reaches  the  waterside  slope,  the  rate  of  lower¬ 
ing  of  the  invert  level  increases  rapidly  and  the  effects 
of  horizontal  flow  contraction,  represented  by  cw2  are 
assumed  to  become  important.  cw  is  now  equal  to  the 
product  of  cw i  and  cw2.  The  third  weir  coefficient,  cw3, 
replaces  the  first  two  once  the  vertical  erosion  process 
has  stopped  and  purely  represents  the  effects  of  the 
horizontal  contraction  at  the  final  breach  stage. 


2.1.1  Surface  erosion 

When  surface  erosion  drives  the  breach  process,  the 
crest  lowers  while  the  landside  slope  retreats  towards 
the  waterside  slope.  The  flow  along  the  landside  slope 
accelerates  and  reaches  maximum  velocity  at  the  bot¬ 
tom  of  the  landside  slope.  With  the  increase  in  flow 
velocity,  the  bed  shear  stresses,  and  corresponding 
erosion  rates,  also  increase.  One  would  expect  this 
to  lead  to  an  increasingly  steepening  slope.  However, 
this  appears  not  to  be  the  case.  It  has  therefore  been 
assumed  that  either  live -bed  scour  effects,  or  geotech¬ 
nical  failures  prevent  the  slope  from  steepening  and 
smooth  out  differences  in  the  erosion  rates  along  the 
landside  slope.  In  AREBA  this  has  been  represented 
by  assuming  that  the  landside  slope  gradient  remains 
equal  to  its  initial  value.  The  overall  erosion  is  as¬ 
sumed  to  be  a  combination  of  the  downward  erosion 
of  the  crest  and  the  retreat  of  the  landside  slope  to¬ 
wards  the  waterside  slope.  The  breach  invert  level 
controls  the  discharge  through  the  breach.  Initially  the 
breach  invert  level  is  solely  influenced  by  the  down¬ 
ward  erosion  of  the  crest.  However  when  the  erosion 
of  the  landside  slope  reaches  the  waterside  slope,  the 
rate  at  which  the  breach  invert  level  drops  increases 
with  increase  of  the  flow  velocities.  The  downward 
and  lateral  erosion  rates  are  quadratically  dependent 
on  the  flow  velocity.  Once  the  water  level  in  the 
breach  falls  below  the  crest  level  the  breach  grows  lat¬ 
erally  as  the  flow  undercuts  the  breach  sides  and  un¬ 
dermines  the  embankment  causing  the  soil  above  the 
undercut  sections  to  become  unstable  and  collapse. 
Most  breach  models  deal  with  the  breach  evolution 
process  by  assuming  a  predefined  breach  shape  which 
is  rectangular,  triangular,  trapezoidal,  or  parabolic. 
The  eventual  breach  geometry  then  follows  from  the 
breach  depth,  and  pre-defined  breach  shape  (Morris 
et  al.,  2009).  The  majority  of  other  models  that  dif¬ 
ferentiate  between  the  breach  erosion  and  breach  side 
slope  stability  processes  deal  with  the  latter  using  a 
simplified  approach.  The  slope  angles  on  both  sides 
of  the  breach  are  usually  found  to  be  steeper  than 
the  corresponding  repose  angle  angle  of  the  embank¬ 
ment  soil.  This  is  because  the  shear  failure  strength  of 
soils  is  related  to  the  amount  of  water  in  the  pores. 
Richards  (1931)  found  a  relationship  between  soil 
pore  pressure  and  moisture  content.  When  a  soil  dries 
out,  the  matric  suction  increases.  As  the  pore  pres¬ 
sure  becomes  progressively  negative  with  a  decrease 
in  moisture  content  the  soil  stability  is  eventually  af¬ 
fected.  The  shear  strength  of  the  soil  depends  on  the 
scale  of  matric  suction.  In  AREBA,  the  effects  of 
the  matric  suction  on  the  side  slope  instability  have 
been  accounted  for  by  assuming  vertical  breach  sides. 
AREBA  requires  a  user-defined  initial  breach  depth 
and  width.  Flow  over  the  intact  landside  slope  is  as¬ 
sumed  to  be  bounded  by  the  same  rectangular  pro- 
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file  as  the  breach.  Hence,  the  flow  is  assumed  not 
to  spread  along  the  landside  slope.  The  breach  width 
is  assumed  as  spatially  constant,  and  the  lateral  and 
downward  erosion  rates  are  assumed  to  depend  solely 
on  the  mean  flow  velocities.  The  lateral  erosion  rate 
has  been  expressed  as  a  function  of  the  downward 
erosion  rate.  For  simplicity  the  foundation  of  the  em¬ 
bankment  is  assumed  to  be  non-erodible.  The  valid¬ 
ity  of  this  assumption  is  questionable  though  erosion 
below  the  bed  level  of  the  reservoir  has  a  negligible 
effect  on  the  discharge  and  hence  the  shape  of  the 
breach  hydrograph.  Once  the  vertical  erosion  reaches 
to  the  foundation  level  and  stops,  the  lateral  erosion 
rate  is  assumed  to  be  a  function  of  the  breach  bed 
shear  stress.  In  AREBA,  the  downward  erosion  rate 
of  the  crest  follows  from 

Ed  =  K  (^pu2UJr  -  ,  (2) 


where  Ed  is  the  downward  erosion  rate  in  m/s,  K  is 
the  soil  erodibility  coefficient  in  m3  /Ns,  p  is  the  den¬ 
sity  of  water  in  kg/m3,  n  the  Manning  coefficient  in 
s/m1'3,  g  the  acceleration  due  to  gravity  in  m/s2,  and 
d  the  water  depth  over  the  crest  in  m  .  The  downward 
crest  erosion  of  the  assumed  trapezoidal  shaped  em¬ 
bankment  leads  the  crest  width  to  increase  with  a  rate 
equal  to  the  downward  erosion  rate  multiplied  with 
factor  S  which  follows  from 


2 bi  tbo 


(3) 


where  ibo  and  ibi  represent  the  waterside  and  land- 
side  slope  gradients.  Recall  the  assumptions  that  the 
landside  slope  remains  straight  and  equal  to  its  ini¬ 
tial  value,  and  that  no  spreading  of  the  flow  occurs 
along  the  landside  slope.  Using  these  assumptions  and 
the  Belanger  equations,  the  shear  stress  has  been  ex¬ 
pressed  as  a  function  of  the  distance  along  the  land- 
side  slope.  With  the  flow  depth  assumed  equal  to  the 
critical  depth  dc  at  the  top  of  the  landside  slope  and 
approaching  the  normal  depth  dn  further  down  the 
landside  slope,  the  erosion  rate  perpendicular  to  the 
slope  Es  becomes 
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where  dn  (■ m )  replaces  the  value  of  the  hydraulic  ra¬ 
dius.  L  is  the  length  in  (  rn)  over  which  a  change  in  the 
normal  depth  Ad  varies  with  a  factor  e,  and  uc  (m/s), 
is  the  velocity  corresponding  to  the  critical  depth  dc 
(m)  that  follows  from 


where  q  is  the  width  averaged  discharge  in  m2  /  s,  and 
g  the  gravitational  constant,  a  follows  from 


a  = 


gl/3n2/3 

dm!J3 


(6) 


The  erosion  perpendicular  to  the  landside  reduces  the 
crest  width,  the  rate  Rc  of  which  follows  from 


(7) 


where  Es  is  the  numerical  averaged  mean  erosion  rate 
along  the  landside  slope,  and  ibi  is  the  landside  slope 
gradient.  From  the  moment  the  landside  slope  reaches 
the  waterside  slope  the  rate  at  which  the  hydraulic 
control  lowers  is  given  by 
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where  ibo  is  the  outward  slope  gradient.  The  cal¬ 
culation  of  the  total  breach  widening  rate  has  been 
adopted  from  HR  BREACH  and  is  assumed  to  be  pro¬ 
portional  to  1.4  times  the  rate  at  which  the  crest  low¬ 
ers  (Mohamed,  2002).  Once  the  crest  height  equals 
the  foundation  level  of  the  embankment  no  further 
vertical  erosion  is  assumed  to  occur.  The  breach 
widening  rate  now  follows  from 
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where  Ew  is  the  rate  of  breach,  or  hydraulic  control, 
widens.  The  breach  stops  growing  when  the  shear 
stresses  exerted  by  the  breach  flow  become  lower  than 
the  critical  shear  stress. 


2.1.2  He adcut  erosion 

Headcut  erosion  is  often  taken  to  be  the  driving  breach 
formation  process  for  an  embankment  composed  of 
cohesive  soil.  Irregularities  along  the  landside  slope 
initiate  a  cascade  formation  that  defines  the  headcut 
erosion  process.  The  cohesive  properties  of  the  soil 
prevent  geotechnical  failures  from  occurring  and  in¬ 
stead  allow  small  cascades  to  initiate  and  grow  in  step 
size.  Where  the  cascade  flow  impacts  on  the  embank¬ 
ment  surface  the  flow  disperses  and  undercuts  the  em¬ 
bankment.  The  local  erosion  rate  depends  upon  the 
energy  with  which  the  flow  hits  the  soil  and  therefore 
on  the  height  of  the  cascades.  At  the  point  when  ero¬ 
sion  causes  the  cascades  to  reach  the  waterside  slope, 
the  flow  rate  increases  rapidly  due  to  the  lowering  of 
the  crest  (i.e.  hydraulic  control)  and  an  increase  in  the 
horizontal  contraction  of  the  flow.  Whilst  the  increase 
in  flow  leads  to  an  increase  in  the  energy  with  which 
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the  flow  hits  the  soil,  the  reduction  in  crest  height  also 
reduces  the  height  over  which  the  flow  drops  before 
hitting  the  soil  and  therefore  reduces  the  flow  energy. 
Zhu  et  al.  (2008)  mention  the  growth  of  headcuts  due 
to  an  increase  in  flow  velocity.  From  field  observa¬ 
tions,  Temple  et  al.  (2005)  assumed  that  the  down¬ 
ward  erosion  of  the  crest  is  negligible  during  the  re¬ 
treat  of  the  headcut  through  the  soil.  In  AREBA,  head- 
cut  failure  is  modelled  using  the  method  described  by 
Temple  et  al.  (2005).  This  method  is  based  on  obser¬ 
vations  of  large  scale  breach  experiments  of  cohesive 
clay  embankments.  The  following  assumptions  made 
by  Temple  et  al.  (2005)  regarding  the  overall  erosion 
process  have  been  adopted  for  the  development  of 
AREBA. 

1 .  The  crest  level  does  not  lower  before  the  headcut 
has  reached  the  upstream  slope. 

2.  The  cross  sectional  breach  shape  is  rectangular. 

3.  The  foundation  of  the  embankment  is  not  erodi- 
ble 

4.  Headcut  initiates  at  the  top  of  the  landside  slope. 

The  first  assumption  is  in  conflict  with  the  physical 
process  as  described  by  Zhu  et  al.  (2008)  who  men¬ 
tion  that  headcuts  grow  due  to  an  increasing  flow  ve¬ 
locity  caused  by  downward  erosion  of  the  crest.  How¬ 
ever  as  mentioned  above,  Temple  et  al.  (2005)  as¬ 
sumed  (based  on  field  observations)  that  the  down¬ 
ward  erosion  of  the  crest  is  negligible  during  the  re¬ 
treat  of  the  headcut  through  the  soil.  This  can  be  ex¬ 
plained  by  the  positive  effect  of  cohesion  on  the  criti¬ 
cal  erosion  shear  stress  of  the  soil.  It  should  be  noted 
that  the  first  assumption  may  limit  the  validity  range 
of  the  model.  For  example,  just  below  the  threshold 
between  headcut  and  surface  erosion  the  critical  shear 
stress  may  be  relatively  low.  Hence  the  assumption 
of  no  vertical  crest  erosion  may  cause  an  error  in  the 
predicted  hydrograph.  The  second  and  third  assump¬ 
tion  have  been  discussed  under  the  header  of  surface 
erosion.  With  regard  to  the  fourth  assumption,  the  rel¬ 
atively  sharp  transition  in  slope  angle  from  the  hori¬ 
zontal  crest  to  the  landside  slope  would  explain  a  sep¬ 
aration  of  the  flow  at  this  stage  that  initiates  the  head- 
cut  erosion  process.  However  the  irregularities,  which 
are  often  present  along  the  landside  slope  could  ini¬ 
tiate  the  headcut  failure  process  anywhere  along  the 
landside  slope  and  the  sharp  transition  at  the  toe  of 
the  slope  often  provides  a  focal  point  for  erosion.  The 
constitutive  equations  for  the  modelling  of  the  head- 
cut  erosion  process  can  be  found  in  the  paper  by  Tem¬ 
ple  et  al.  (2005)  and  will  not  be  further  addressed  here. 


2.2  Piping 

Piping  is  defined  here  as  the  internal  erosion  of  em¬ 
bankment  material  by  seepage  flow.  Sellmeijer  (1988) 
reports  that  piping  may  stop  once  some  erosion  has 
taken  place  due  to  the  larger  permeability  of  the  em¬ 
bankment  where  the  soil  has  eroded,  compared  with 
the  permeability  of  the  surrounding  soil.  The  increase 
in  permeability  reduces  the  hydraulic  gradient  that 
initiates  the  piping  process.  Piping  failure  is  modelled 
using  the  method  proposed  by  Mohamed  (2002)  and 
applied  in  the  HR  BREACH  model  whereby  an  ini¬ 
tial  open  pipe  connection  has  been  assumed  between 
both  sides  of  the  embankment.  The  erosion  along  the 
perimeter  of  the  pipe  causes  the  pipe  diameter  to  in¬ 
crease  to  the  point  that  the  soil  above  the  pipe  fails 
and  slumps  downs  into  the  pipe.  Depending  upon  the 
flow  velocities  through  the  pipe  and  the  volume  of 
soil  slumping  down,  the  slumped  material  could  ei¬ 
ther  be  transported  almost  instantaneously  or  could 
cause  the  pipe  flow  to  stop  and  the  embankment  crest 
to  lower  after  which  an  overflow  failure  could  initiate. 
In  AREBA  the  latter  of  the  two  is  assumed  to  occur. 
Soil  properties  like  cohesion,  matric  suction  and  inter¬ 
nal  friction  angle  determine  when  the  material  above 
the  pipe  fails.  Positive  arching  of  the  soil  above  the 
pipe,  whereby  part  of  the  load  above  the  pipe  is  trans¬ 
ferred  to  other  soil  prisms,  reduces  the  load  on  the 
pipe.  The  effects  of  soil  parameters  are  difficult  to  pre¬ 
dict  using  rapid  methods.  For  that  reason  AREBA  as¬ 
sumes  that  soil  cohesion  and  arching  solely  affect  the 
instant  at  which  slumping  occurs.  It  should  be  noted 
that  though  this  is  quite  a  coarse  assumption,  its  ef¬ 
fects  are  expected  to  lie  within  the  bounds  of  uncer¬ 
tainly  associated  with  the  soil  parameters  themselves 
since  the  properties  of  soils  in  embankments  can  be 
quite  uncertain.  The  pipe  is  assumed  to  grow  at  an 
equal  rate  along  its  perimeter  and  to  maintain  its  cir¬ 
cular  cross-sectional  shape.  The  pipe  cannot  grow  be¬ 
low  the  foundation  level  of  the  embankment  or  above 
the  water  level  upstream,  in  which  the  shape  takes  the 
form  of  a  segment.  The  soil  immediately  above  the 
pipe  is  assumed  to  fail  once  its  weight  exceeds  the 
stabilising  forces  along  the  sides  of  the  soil  mass.  The 
steady  pipe  flow  rate  is  calculated  using  Bernouilli’s 
energy  equation  (Nortier  and  De  Koning,  1991) 


Qb  —  A 


2gAH 

hf 


(10) 


where  A  is  the  wet  cross  sectional  area  of  the  pipe  in 
m2,  AH  is  the  head  difference  over  the  pipe  in  m, 
and  hf  is  a  head  loss  factor  representing  contraction 
and  friction  losses.  In  AREBA  the  head  loss  factor  is 
given  by 

hf  —  1.05  +  (11) 

Up 
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where  L  is  the  pipe  length  in  m,  and  Dp  is  the  pipe 
diameter  in  m.  f  is  a  friction  coefficient  determined 
from 

/  D  \ 

f  =  0.2165  (-^°J  (12) 

where  D50  ,  is  the  grain  diameter  in  m  exceeded  by 
50%  of  the  soil  mass.  Due  to  the  long  straight  pipe 
geometry,  the  energy  losses  are  dominated  by  the  wall 
friction  and  hence  the  losses  due  to  outflow  from  the 
pipe  are  neglected.  The  cross  sectional  area  of  the 
pipe  A  is  calculated  assuming  a  circular  shaped  pipe 
that  does  not  grow  beyond  the  embankment’s  founda¬ 
tion  level.  Once  the  pipe  bottom  reaches  the  embank¬ 
ment’s  foundation,  the  pipe  shape  becomes  that  of  a 
truncated  circle.  The  surface  area  of  the  pipe  then  fol¬ 
lows  from 


A  =  -^7r -  ^_D2  +  h2ptana  (13) 


in  which 


a  = 


arccos  (2  hp) 
Dp 


(14) 


where  hp  is  the  distance  between  the  centreline  of  the 
pipe  and  the  foundation  level  of  the  embankment.  The 
pipe  is  assumed  to  grow  uniformly  along  its  length 
according  to 


E  =  2 K  ( phRiw  —  rc)  (15) 


where  R  is  the  hydraulic  radius  of  the  pipe,  which  for 
a  submerged  circular  pipe  is  equal  to  Dp/A ,  with  Dp 
representing  the  pipe  diameter.  iw  is  the  gradient  in 
water  level  over  the  pipe  and  follows  from 
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In  AREBA  the  assumption  has  been  adopted  from 
Mohamed  (2002)  that  the  soil  mass,  above  the  pipe 
fails  when  the  weight  of  the  soil  above  the  pipe  ex¬ 
ceeds  the  friction  at  its  sides.  Mohamed  gives  the  fol¬ 
lowing  relationship  for  the  factor  of  safety  ( FOS ) 
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cAs 

W 


(17) 


where  c  is  the  cohesion  of  the  soil  in  N/rn2,  As  is  the 
cross-sectional  area  of  the  embankment  in  m2  above 
the  centreline  of  the  pipe.  For  simplicity  the  weight 
of  the  soil  W  in  N  is  purely  determined  for  the  soil 
above  the  top  of  the  pipe  from. 


W  =  fapgDpAs  (18) 

where  As  is  the  surface  area  of  the  sides  of  the 
soil  wedge  calculated  from  the  top  of  the  pipe  up¬ 
wards,  Mohamed  gives  the  following  relationship  for 


fa  which  is  an  empirically  based  reduction  factor  rep¬ 
resenting  the  arching  of  the  soil. 

fa  =  9.860  x  10_5r2  -  2.918  x  10"4r  +  3.757  x  10“4 

(19) 

where  r  represents  the  ratio  of  the  distance  between 
the  top  of  the  pipe  and  the  crest,  and  the  pipe  diam¬ 
eter.  The  ratio  r  has  a  minimum  value  of  0.1  and  a 
maximum  value  of  1.6.  For  FOS  <  1,  the  pipe  col¬ 
lapses,  and  the  embankment  crest  is  assumed  to  lower 
a  distance  equal  to  the  height  needed  to  fill  the  vol¬ 
ume  of  the  pipe.  From  this  time  onwards  the  embank¬ 
ment’s  failure  mechanism  transfers  to  that  of  a  surface 
erosion  process. 

3  VERIFICATION  AND  VAFIDATION 
The  results  from  AREBA  presented  and  discussed  in 
this  Section  were  obtained  with  model  run  speeds  of 
less  than  one  second  per  run.  This  makes  AREBA 
suitable  for  use  within  system  risk  models  where 
many  hundreds  or  thousands  of  simulations  may  be 
necessary  to  calculate  the  overall  system  flood  risk. 
As  a  benchmark  for  the  modelling  accuracy,  AREBA 
has  been  tested  against  the  HR  BREACH  model, 
which,  together  with  the  SIMBA  model,  was  selected 
by  the  Dam  Safety  Interest  Group  (DSIG)  breach 
modelling  project  as  models  with  the  potential  for 
further  development  for  industry  use  (Wahl,  2007). 
HR  BREACH  version  4.1  contains  the  headcut  failure 
mode  from  SIMBA,  and  so  is  used  for  benchmarking. 

Validation  of  HR  BREACH  against  different 
datasets  led  to  the  incorporation  of  an  empirical  co¬ 
efficient  of  1.6  for  the  downward  erosion  rate.  This 
factor  was  based  on  model  performance  for  a  range 
of  material  and  structure  types.  To  enable  proper  com¬ 
parison  between  AREBA  and  HR  BREACH  the  same 
coefficient  has  been  adopted  in  AREBA.  Table  1  lists 
the  input  parameters  that  were  kept  constant  for  all 
verification  runs  with  HR  BREACH,  and  Table  2 
gives  the  input  parameters  that  were  varied.  hcrest 
refers  to  the  crest  height,  and  Slope  refers  to  both 
the  landside  as  the  waterside  slope  gradient.  A  box 
shaped  reservoir  with  continuous  constant  inflow  was 
chosen  for  the  upstream  boundary  condition.  An  in¬ 
finitely  large  floodplain  was  chosen  for  the  down¬ 
stream  boundary  condition  to  facilitate  a  clear  com¬ 
parison  between  the  two  models.  In  AREBA’s  surface 
erosion  mode  the  landside  slope  gradient  is  assumed 
to  remain  equal  to  its  initial  value  when  the  land- 
side  slope  retreats  towards  the  waterside  slope.  In  HR 
BREACH  unrealistic  landside  slope  gradients  are  pre¬ 
vented  from  occurring  by  averaging  the  erosion  rates 
at  the  top  of  the  landside  slope  and  downstream  side 
of  the  crest,  and  averaging  the  erosion  rate  at  the  toe 
of  the  landside  slope  and  the  foundation  level.  Fig¬ 
ures  2,  and  3  show  respectively  the  results  of  Run  1 
and  Run  2.  whereby  the  steepness  of  the  slope  was 
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Table  1 :  Fixed  input  parameters  HR  BREACH  validation 
runs _ 


Parameter 

Value 

Unit 

Initial  water  level  reservoir 

hcrest  ~  0.3  * 

m 

Crest  width 

6 

m 

Discharge  into  reservoir 

1 

m3  /  s 

Surface  area  reservoir 

30000 

2 

rn 

Initial  breach  depth 

0.25 

m 

Crest  length 

200 

rn 

Soil  cohesion 

1  x  10"5 

kN/m 2 

Soil  tension  strength 

1  x  10"5 

kN/m 2 

Friction  angle 

30 

degrees 

Porosity 

0.2 

- 

Dry  unit  weight 

16.8 

kN/m3 

Critical  shear  stress 

1  x  10"5 

N/m 2 

Manning  coefficient 

0.035 

s/m3 

Weir  coefficient 

1 

- 

*  hcrest  =  crest  height 


Table  2:  Variable  input  parameters  HR  BREACH  valida- 
tion  runs _ 


Run  no. 

Erodibility  coefficient 

(m3  /Ns) 

crest  height 

(m) 

Slope 

m/m 

1 

5  x  10"b 

12 

1/3 

2 

5  x  10-6 

6 

1/3 

3 

5  x  10-6 

9 

1/6 

4 

5  x  10-6 

9 

1/2 

kept  constant.  Differences  in  approaches  to  limit  the 


Figure  2:  Hydrograph  output  HR  BREACH  and  AREBA: 
Run  1 


steepness  of  the  landside  slope  lead  to  differences  in 
model  outcomes,  which  are  evident  by  comparing  the 
results  from  Runs  3  and  4  which  only  differ  in  the 
steepness  of  the  landside  slope  (see  Figures  4  and  5). 
To  get  a  better  estimate  of  the  impact  of  the  assump¬ 
tions,  AREBA  has  been  validated  against  the  EU  IM¬ 
PACT  project  experiments  for  the  surface  erosion  fail- 


Figure  3:  Hydrograph  output  HR  BREACH  and  AREBA: 
Run  2 


ure  and  piping  failure  (IMPACT,  2005). 

3.1  Validation 

For  proper  validation  against  the  EU  IMPACT  exper¬ 
iments  AREBA  was  adapted  to  use  a  hydrograph,  and 
stage-area  curve  as  boundary  conditions.  Hassan  and 
Morris  studied  the  accuracy  of  the  input  and  output 
values  of  the  IMPACT  project  experiments  (Hassan 
and  Morris,  2009).  Besides  being  uncertain  about  the 
exact  values  to  use  as  input,  some  of  the  model  in¬ 
put  parameters  were  undefined  and  needed  to  be  es¬ 
timated.  The  problem  with  the  uncertainty  of  the  in¬ 
put  parameters  is  that  more  unknowns  exist  than  there 
are  validation  tests  available.  Several  combinations 
of  input  parameters  can  therefore  be  found  that  all 
give  a  near  perfect  match  with  the  validation  data. 
This  does  not  however  mean  that  the  model  repre¬ 
sents  the  physical  processes  correctly.  Therefore  in¬ 
stead  of  calibrating  the  unknown  model  parameters 
to  fit  the  experimental  data,  upper  and  lower  limits 
have  been  estimated  for  the  input  parameters  based  on 
the  data,  reports,  and  empirical  relationships.  These 
limits  formed  the  boundaries  of  the  uniform  distribu¬ 
tion  from  which  input  values  were  drawn.  By  running 
AREBA  500  times  for  different  combinations  of  input 
parameters  an  envelope  was  obtained  between  which 
the  experimental  results  are  expected  to  lie.  After  hav¬ 
ing  obtained  the  outer  envelop,  the  same  process  was 
repeated  for  half  the  range  of  input  values  while  main¬ 
taining  the  same  mean  value  in  order  to  study  the 
impact  on  the  bounds  of  not  knowing  the  input  val¬ 
ues.  Table  3  shows  the  upper  and  lower  limits  and  the 
mean  values  used  for  validation  of  AREBA’s  surface 
erosion  mode,  where  n  is  the  Manning  coefficient,  rc 
is  the  critical  shear  stress,  K  is  the  erodibility  of  the 
embankment  soil,  cw i  and  cw2  are  respectively  the  first 
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Figure  4:  Hydrograph  output  HR  BREACH  and  AREBA:  Figure  5:  Hydrograph  output  HR  BREACH  and  AREBA: 

Run  3  Run  4 


Table  3:  Range  of  input  parameters  AREBA  for  IMPACT 
surface  erosion  validation _ 


Symbol 

Lower  bound 

Mean 

Upper  bound 

n 

0.025 

0.035 

0.045 

Tc 

1.8 

4.2 

7.1 

K 

1.8  x  10-7 

9  x  10-6 

1.8  x  10"5 

Cwl 

0.9 

1 

1.1 

Cw2 

1 

1.2 

1.4 

sw 

1/1.7 

1/1.8 

1/1.9 

Si 

1/1.6 

1/1.65 

1/1.7 

hbr 

0.1 

0.105 

0.12 

and  second  discharge  coefficient,  Sw  and  Si  are  re¬ 
spectively  the  waterside  and  landside  slope  gradient, 
and  hbr  the  initial  breach  depth.  The  third  weir  coeffi¬ 
cient  cw3  was  set  to  be  1.  Figure  6  displays  the  maxi¬ 
mum  and  minimum  values  of  the  discharge  found  per 
run.  Figure  7  shows  the  breach  width  development. 

The  term  large  range  mentioned  in  the  figures  rep¬ 
resents  the  results  found  by  applying  the  Lower  and 
Upper  bound  as  bounds  for  the  uniform  distribution. 
The  term  small  range  refers  to  the  outcome  obtained 
when  using  the  average  of  the  mean  and  the  lower 
bound,  and  the  average  of  the  mean  and  upper  bound 
as  respectively  the  lower  and  upper  bound  values.  The 
final  breach  width  measured  during  the  IMPACT  ex¬ 
periment  was  10m  which  falls  between  the  upper  en¬ 
velopes  given  by  AREBA  for  the  large  range  and  the 
small  range.  The  lower  envelopes  found  with  AREBA 
coincide  with  the  time-axes  in  Figures  6  and  7  and 
therefore  have  not  been  plotted. 

The  technical  report  describing  the  IMPACT  pip¬ 
ing  experiment  mentions  a  rapid  pipe  failure  after  ap¬ 
proximately  20  minutes.  The  moment  of  opening  the 
pipe  was  estimated  from  the  field  data  and  was  ex¬ 
pected  to  coincide  with  a  sudden  drop  in  water  level 


(see  Figure  8).  The  exact  moment  at  which  the  pipe 
was  opened,  relative  to  the  start  time  of  the  given 
inflow  hydrograph,  is  however  unknown.  Before  the 
drop  in  water  level,  the  reservoir  was  filled  slowly. 
During  the  experiment  it  was  attempted  to  keep  the 
upstream  water  level  constant  irrespective  of  the  out¬ 
flow.  The  ability  to  control  the  upstream  water  level 
by  accounting  for  the  effects  of  the  outflow  from  the 
pipe  had  strongly  been  underestimated  and  caused  for 
a  rapid  increase  of  inflow  into  the  reservoir.  Validation 
of  AREBA  against  the  IMPACT  field  data  was  com¬ 
plicated  due  to  the  unknown  opening  time  of  the  pipe 
relative  to  the  given  inflow  hydrograph,  the  unknown 
discharge  into  the  reservoir  before  opening  the  pipe, 
and  the  sudden  increase  in  discharge  into  the  reser¬ 
voir  just  after  opening  the  pipe.  Treating  the  relative 
opening  time  of  the  pipe  as  a  statistical  variable  in 
the  analysis  resulted  in  a  very  large  uncertainty  bound 
since  the  sudden  increase  in  discharge  just  after  open¬ 
ing  the  pipe  resulted  in  widely  varying  outcomes  for 
different  initiation  times.  A  similar  problem  was  en¬ 
countered  with  specifying  a  wide  range  for  the  credi¬ 
bility.  The  range  for  the  credibility  was  therefore  lim¬ 
ited.  The  inflow  before  the  pipe  was  opened,  and  the 
time  at  which  the  pipe  opened  were  obtained  by  fit¬ 
ting  the  data  slightly  to  the  found  outcome  and  were 
respectively  estimated  at  0.6m3/s  and  13450s.  The 
range  of  each  of  the  input  parameters  is  given  in  Ta¬ 
ble  4.  The  model  outcomes  are  presented  by  Figures 
9  and  8.  Due  to  the  small  statistical  range  of  the  vari¬ 
ables  the  piping  case  was  only  run  for  the  wide  range 
in  data. 

4  DISCUSSION 

The  accuracy  with  which  AREBA  predicts  a  flood  hy¬ 
drograph  is  directly  related  to  how  well  the  underlying 
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Figure  6:  Maximum  and  minimum  discharge  (per  time 
step)  predicted  by  AREBA  for  500  runs,  compared  against 
the  IMPACT  field  data 


Table  4:  Range  of  input  parameters  AREBA  for  IMPACT 
piping  validation  case 


Symbol 

Lower  bound 

Mean 

Upper  bound 

n 

0.025 

0.035 

0.045 

c 

18 

20 

22 

P 

2000 

2100 

2200 

rc 

0 

5 

10 

K 

9e  —  6 

le  —  5 

l.le -5 

t'w  1 

0.9 

1 

1.1 

Cw2 

1 

1.2 

1.4 

sw 

1/1.3 

1/1.35 

1/1.4 

Si 

1/1.3 

1/1.35 

1/1.4 

hp 

0.29 

0.3 

0.31 

Merest 

4.3 

4.4 

4.5 

WCrest* 

2.8 

2.9 

3 

*  Wcrest  =  crest  width 


assumptions  represent  the  physics  of  the  breach  for¬ 
mation.  The  impact  of  the  error  made  assuming  one¬ 
dimensional  flow  behaviour,  for  example,  is  quadrati- 
cally  amplified  when  calculating  the  erosion  rate  due 
to  the  quadratic  dependence  of  the  shear  stress  on  the 
mean  flow  velocity.  The  assumptions  that  influence 
the  accuracy  of  the  hydrograph  prediction  are 

•  Negligible  flow  accelerations  perpendicular  to 
the  main  flow  direction 

•  Horizontal  flow  contraction  commences  only 
when  the  landside  slope  retreat  has  reached  the 
waterside  slope 

•  Flow  velocities  in  the  reservoir  are  negligible 

•  The  breach  width  is  spatially  constant. 


Figure  7 :  Maximum  and  minimum  breach  width  (per  time 
step)  predicted  by  AREBA  for  500  runs  compared  against 
the  IMPACT  field  data 

•  The  flow  depth  on  top  of  the  landside  slope  is 
equal  to  the  critical  depth 

•  The  bed  roughness  is  spatially  constant 

•  The  critical  shear  stress  is  spatially  constant 

These  assumptions  originate  from  the  objective  to 
have  a  high  computational  speed  but  limit  the  abil¬ 
ity  of  AREBA  to  simulate  the  physics  accurately.  The 
physics  behind  breaching  is  in  itself  still  relatively 
inadequatly  understood.  The  question  as  to  how  to 
change  the  assumptions  for  a  more  accurate  descrip¬ 
tion  of  the  breach  formation,  while  maintaining  a  fast 
model,  is  difficult  to  answer.  For  example,  the  Man¬ 
ning  coefficient  is  known  to  vary  along  the  breach,  but 
it  is  unknown  where  and  by  how  much  it  varies,  and 
is  consequently  assumed  as  constant. 

Verification  of  AREBA  against  HR  BREACH  indi¬ 
cated  a  very  close  agreement  between  the  models  for 
Run  2  (Figure  3).  The  input  of  Run  2  differs  from  that 
of  Run  1  (Figure  (2)  in  the  ratio  of  the  crest  height 
over  the  crest  width.  In  AREBA,  the  crest  is  assumed 
to  erode  downwards  while  the  landside  slope  retreats 
towards  the  waterside  slope.  For  relatively  high  ra¬ 
tios  of  the  crest  width  to  the  crest  height,  or  rela¬ 
tively  high  values  for  the  erodibility,  the  crest  will  al¬ 
ready  have  eroded  down  to  the  foundation  level  be¬ 
fore  the  retreat  of  the  landside  slope  reaches  the  wa¬ 
terside  slope.  Figure  5  shows  a  nice  fit  is  obtained 
between  both  model  outcomes  for  a  landside  slope 
of  1:2.  To  the  contrary  Figure  4  shows  less  agree¬ 
ment  for  a  landside  slope  of  1:6.  It  should  be  noted 
that  the  model  outcomes  of  AREBA  are  identical  for 
the  two  landside  slope  gradients  because  the  down¬ 
ward  crest  erosion  has  already  reached  the  embank¬ 
ment’s  foundation  before  the  landside  slope  retreat 
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Figure  8:  Maximum  and  minimum  breach  width  (per  time 
step)  predicted  by  AREBA  for  500  runs  compared  against 
the  IMPACT  field  data 


has  reached  the  waterside  slope.  Despite  the  fact  that 
AREBA  and  HR  BREACH  disagree  in  the  value  for 
the  peak  discharge,  they  do  agree  on  the  overall  shape 
of  the  hydrograph  and  the  time  of  the  peak  discharge. 
Both  assumptions  for  the  prevention  of  unrealistic 
steep  slopes  are  empirical.  Without  modelling  all  the 
physics  involved  it  is  not  possible  to  state  which  as¬ 
sumption  has  the  greater  effect  on  the  accuracy  to 
which  the  physical  process  is  represented. 

To  validate  AREBA,  the  influence  of  the  uncer¬ 
tainties  in  the  input  parameters  was  assumed  to  out¬ 
weigh  the  effects  of  simplifying  the  physics  in  the 
constitutive  equations.  AREBA  was  validated  against 
field  data  using  a  range  of  input  parameters.  Figures 
6  and  7  show  that  for  the  surface  erosion  case,  the 
hydrograph  and  final  breach  widths  obtained  from 
field  measurements  lay  within  the  envelopes  found  by 
AREBA  for  the  large  range,  and  close  to  envelopes 
found  for  the  small  range.  Since  the  model  is  volume 
conservative  the  model  is  capable  of  reproducing  the 
outcome  of  the  field  data.  However,  it  is  unknown  if 
the  close  match  is  only  found  by  using  values  for  the 
Manning  coefficient  and  erodibility  coefficient  that 
do  not  represent  the  actual  situation,  and  hence  it  is 
unknown  to  what  accuracy  AREBA  is  able  to  sim¬ 
ulate  the  physical  breach  processes.  Similar  to  other 
breach  models  the  low  quality  of  the  limited  amount 
of  datasets  available  makes  it  difficult  to  make  a  state¬ 
ment  regarding  the  capability  of  the  AREBA  model 
to  simulate  the  physical  processes  of  an  embankment 
failure  correctly. 

The  nature  of  the  EU  IMPACT  piping  experiment 
did  not  permit  validation  based  on  a  large  range  of 
input  parameters.  The  sudden  increase  in  discharge 


Figure  9:  Maximum  and  minimum  discharge  (per  time 
step)  predicted  by  AREBA  for  500  runs,  compared  against 
the  IMPACT  field  data 

into  the  reservoir  combined  with  an  unknown  time  at 
which  the  pipe  was  opened  resulted  in  model  condi¬ 
tions  that  caused  the  model  to  be  highly  sensitive  for 
the  input  parameters.  For  example,  an  under-predicted 
Manning  coefficient  or  erosion  coefficient  leads  to  a 
smaller  growth  in  pipe  diameter  and  a  smaller  out¬ 
flow.  A  sudden  increase  in  inflow  in  the  reservoir  then 
results  in  unrealistic  high  water  levels  in  the  reservoir, 
which  at  the  instant  the  pipe  would  fail  would  cause 
for  unrealistic  outcomes  when  the  model  changes  to 
surface  erosion  mode.  Both  Figures  9  and  8  show  a 
second  spike  in  the  model  outcomes  which  results 
from  the  assumption  that  when  the  pipe  fails  the  body 
of  soil  above  the  pipe  fills  up  the  volume  of  the  pipe. 
During  the  experiment  the  soil  that  slumped  in  the 
pipe  was  washed  away  instantly.  However  since  the 
instantaneous  removal  of  slumped  material  is  case  de¬ 
pendent,  the  description  of  the  physics  by  the  model 
has  not  been  changed. 


As  shown  by  the  dataset  used  for  validation  the  soil 
parameters  are  often  poorly  known  or  unknown.  In 
some  cases  even  the  type  of  construction  material  is 
unknown.  The  contribution  of  modelling  the  influence 
of  soil  parameters  on  the  breach  process  of  an  em¬ 
bankment  is  therefore  variable.  Even  if  the  soil  pa¬ 
rameters  were  known,  it  is  still  poorly  understood  as 
to  how  these  interact  with  the  hydrological  aspects. 
The  assumptions  used  in  AREBA  are  crude  and  do  af¬ 
fect  the  ability  of  the  model  to  simulate  a  breach  cor¬ 
rectly  though  the  limited  amount  of  input  data  avail¬ 
able  make  the  use  of  highly  sophisticated  models  over 
the  simplified  ones  questionable. 
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5  CONCLUSIONS 

AREBA  is  able  to  rapidly  predict  a  flood  hydrograph, 
breach  width,  and  breach  depth  development  for  sur¬ 
face,  headcut,  and  piping  failures  for  trapezoidal 
shaped  homogeneous  embankments.  The  speed  at 
which  the  AREBA  model  runs  (less  than  Is)  makes  it 
suitable  for  an  application  in  system  risk  models.  The 
time  of  occurrence  of  the  peak  discharge  corresponds 
closely  with  those  found  by  HR  BREACH.  Any  dif¬ 
ference  in  the  prediction  of  the  peak  discharge  was 
found  to  originate  from  assumptions  with  respect  to 
preventing  unrealistic  values  for  the  steepness  of  the 
landside  slope.  However,  in  order  to  better  assess  the 
accuracy,  and  validity  of  the  underlying  assumptions 
of  AREBA  and  other  breach  models,  more  accurate, 
large  scale  experiments  are  required.  At  this  point  it  is 
possible  to  state,  based  on  the  validation  of  AREBA 
against  the  IMPACT  experiments  that  AREBA  is  able 
to  reproduce  the  shape  of  the  breach  hydrograph  for 
input  values  that  lie  within  the  range  of  uncertainty. 
In  comparison  to  the  use  of  existing  simple  breach  as¬ 
sumptions  within  system  risk  models,  or  even  the  use 
of  regression  equations  that  predict  the  potential  peak 
discharge  from  a  breach,  the  AREBA  model  should 
provide  a  more  accurate  answer  and  in  doing  so,  also 
provides  an  estimate  of  the  breach  flood  hydrograph 
shape. 
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Abstract 

Effective  management  of flooding  requires  models  that  are  capable  of  quantifying 
flood  risk.  Quantification  of food  risk  involves  both  the  quantification  of 
probabilities  of flooding  and  the  associated  consequences.  Modern  flood  risk  models 
account  for  the  probabilities  of  extreme  hydraulic  loading  events  and  also  include  a 
probabilistic  representation  of  the  performance  of flood  defence  infrastructure  and  its 
associated  reliability.  The  spatial  and  temporal  variability  of  flood  events  makes 
probabilistic  representation  of  the  hydraulic  loading  conditions  on  the  flood  defences 
complex.  In  the  system  method  used  widely  within  England  and  Wales,  simplifying 
assumptions  relating  to  the  spatial  dependence  of  flood  events  are  made.  Recent 
research  has  shown  the  benefits  of  using  improved  multivariate  extreme  value 
methods  to  define  the  hydraulic  loading  conditions  for  flood  risk  assessment.  This 
paper  describes  the  development  of  a  new  modelling  system  that  improves  the  systems 
based  risk  analysis  model  currently  applied  in  practice,  through  the  incorporation  of 
a  multivariate  extreme  value  model.  The  new  system  has  been  presented  on  a  case 
study  site  in  the  North  West  of  England. 

Keywords 

Flood  risk  analysis,  Spatial  dependence,  Multivariate  extremes,  Systems  modelling, 
Flood  defence  failure,  Reliability  analysis 


Introduction 

Flooding  is  a  global  problem  with  extreme  flood  events  having  major  consequences 
recorded  at  different  locations  relatively  frequently.  It  has  long  been  recognised  that 
risk  management  approaches  offer  a  number  of  benefits  over  more  traditional 
detenninistic  design  event  based  approaches,  USAGE  (1996),  Sayers  et  al.  (2002),  for 
example.  A  primary  component  of  successful  flood  risk  management  is  the  ability  to 
quantify  flood  risk  and  to  simulate  the  risk  reduction  that  results  from  introducing 
different  mitigation  measures. 

Flood  risk  is  generally  regarded  as  a  function  of  probability  and  consequence  and 
quantifying  flood  risk  can  be  complex.  Extreme  flood  events  at  any  particular  spatial 
location  occur,  by  definition,  rarely  and  there  is  limited  data  with  which  to  verify 
models.  The  temporal  variability  of  flooding  can  also  be  complex.  Some  flash  flood 
events  last  for  a  short  period  of  time,  a  day,  for  example,  whilst  others  can  last  for 
many  weeks.  Flooding  can  be  localised  in  terms  of  its  impact  or  can  effect  large 
spatial  scales.  Within  the  UK  for  example,  the  flooding  of  the  village  of  Boscastle 
(August,  2004),  that  took  place  over  a  day,  Roca-Collel  and  Davison  (2010),  can  be 
contrasted  with  the  summer  floods  of  2007,  Marsh  and  Hannaford  (2007),  that  lasted 
for  a  period  of  six  weeks  and  affected  most  of  England  and  parts  of  Wales.  Flooding 
can  arise  from  different  sources;  pluvial,  fluvial,  waves  and  surges  for  example. 

These  different  sources  can  often  interact  to  exacerbate  the  impact  of  the  flood  event. 

As  risk  is  a  function  of  probability,  establishing  the  joint  probability  of  different 
combinations  of  flood  inducing  events  arising  has  been  the  subject  of  much  research. 
Examples  include  extreme,  fluvial  flows  at  confluences,  estuarine  water  levels, 
astronomical  tides  and  metrological  surges,  waves  and  sea  levels,  see  for  example, 
Tawn  and  Vassie  (1989),  Acreman  (1994),  Bruun  and  Tawn  (1998),  Hawkes  et  al. 
(2002),  Cai  et  al.  (2008),  Wahl  et  al.  (2012).  These  approaches  have,  in  the  past, 
tended  to  have  been  restricted  in  tenns  of  the  number  of  variables  that  have  been 
considered  and  the  spatial  extents  that  are  covered.  The  reasons  for  this  primarily 
relate  to  constraints  placed  by  the  statistical  models  that  have  been  used.  In  particular, 
traditional  multivariate  extreme  value  statistical  models  have  constraints  relating  to 
the  handling  of  the  dependence  structure.  These  constraints  make  extensions  to  large 
spatial  scales,  where  the  degree  of  dependence  between  extremes  of  a  given  variable 
can  cover  a  wide  range,  impractical  to  implement.  Recent  developments  in 
multivariate  extreme  value  methods,  Heffernan  and  Tawn  (2004),  have  however, 
removed  some  of  these  constraints  and  opened  opportunities  for  improving  flood  risk 
analysis  methods. 

Risk  analysis  models  of  flooding  systems  that  incorporate  the  likelihood  of  flood 
defence  failure  are  now  recognised  as  an  effective  means  of  supporting  the  flood  risk 
management  process.  These  models  have  been  applied  at  national  and  regional  scales 
for  supporting  decisions  relating  to  long  tenn  climate  change  adaptation,  strategic 
planning  and  asset  management  prioritisation  purposes,  see  US  ACE  (1996), Hall  et  al. 
(2006)  ,  Evans  et  al.  (2006),  Apel  et  al.  (2004),  Gouldby  et  al.  (2008b)  and 
Vorogushyn  et  al.  (2010),  Woodward  et  al.  (2011),  Woodward  et  al.  (2012),  for 
example.  Within  England  and  Wales,  the  Environment  Agency  has  been  utilising 
these  flood  defence  systems  based  models  for  around  a  decade,  Environment  Agency 


(2003),  Hall  et  al.  (2003),  Gouldby  et  al.  (2008a).  The  latter  methodology  has  now 
been  implemented  within  the  Environment  Agency's  flood  risk  modelling  decision 
support  system  (MDSF2),  Environment  Agency  (2011a),  where  it  is  being  routinely 
applied  for  systems  based  flood  risk  analysis.  One  of  the  limitations  of  this  current 
method  however,  relates  to  the  handling  of  dependencies  of  flood  events  over  large 
spatial  scales.  To  address  this  problem  the  Environment  Agency  commissioned  a 
research  project  to  further  explore  aspects  of  spatial  dependence  in  relation  to  flood 
risk  analysis.  More  specifically,  the  application  of  the  multivariate  extreme  value 
statistical  model  of  Heffeman  and  Tawn  (2004)  was  explored  in  the  context  of  flood 
risk  analysis.  The  findings  of  this  research  project  have  been  published  in  a  series  of 
reports,  Environment  Agency  (2011b),  Environment  Agency  (2011c),  Environment 
Agency  (201  Id),  and  also  by  Lamb  et  al.  (2010). 

Lamb  et  al.  (2010)  and  Environment  Agency  (2011b)  note  the  potential  for 
integrating  the  statistical  method  of  Heffeman  and  Tawn  (2004)  with  the  system- 
based  risk  model  used  by  the  Environment  Agency,  and  provide  a  thorough 
discussion  of  the  issues.  This  paper  forms  the  natural  extension  of  this  discussion. 
The  Heffeman  and  Tawn  (2004),  multivariate  extreme  value  method,  refined  by  Keef 
et  al.  (2009),  has  been  used  to  generate  boundary  conditions  for  the  system  risk 
model,  described  by  Gouldby  et  al.  (2008a).  The  new  modelling  system  has  been 
applied  on  a  case  study  in  the  North  West  of  England,  with  the  results  described 
below. 


Background  to  system  flood  risk  analysis 

System  risk  models  that  are  currently  applied  in  practice  typically  define  risk  through 
consideration  of  the  aleatory  uncertainty  associated  with  the  random  nature  of 
extreme  flood  events  and  the  epistemic  uncertainty  associated  with  the  potential  for 
flood  defence  infrastructure  to  suffer  failure.  Whilst  there  are  many  other  sources  of 
uncertainty,  Hall  and  Solomatine  (2008),  approaches  to  quantify  these,  Merz  and 
Thieken  (2009b),  Gouldby  et  al.  (2010),  for  example,  are  not  commonly  applied  in 
practice,  primarily  due  to  the  computational  burden  associated  with  the 
implementation  of  the  methods.  The  primary  components  of  the  risk  models  are: 
hydraulic  loads,  described  by  extreme  value  distributions,  flood  defence  infrastructure 
(dykes,  levees  or  embankments  for  example),  defined  by  fragility  curves,  flood 
inundation  simulation  representation  and  functions  that  relate  the  simulated  floods  to 
consequences.  The  flood  defences  are  typically  defined  as  discrete  lengths  with  a 
specific  fragility  curve  prescribed  for  each  length.  The  models  that  are  applied  in 
practice  within  England  and  Wales  assume  the  perfonnance  of  each  of  the  defence 
lengths  is  independent  and  hence  the  risk,  expressed  in  terms  of  the  Expected  Annual 
Damage  (EAD),  is  given  by: 


EAD  =  jfjP(di\X)fx(X)g(di,X)dX  (1) 

z=l 

where  n  is  the  number  of  defence  lengths, ^  is  the  joint  density  of  hydraulic  loads  X 
over  the  defence  lengths,  dt  is  the  defence  system  state  (a  vector  of  length  n  that 
comprises  a  representation  of  the  state  as  failed  or  undamaged  of  each  defence  in  the 
system)  and  g  is  a  function  that  relates  the  defence  system  state  and  hydraulic  loads  to 
the  economic  damage  to  property.  The  sum  in  i  covers  all  2"  possible  configurations 
of  the  defence  system  state  d,. 


The  derivation  of  the  joint  density  of  the  hydraulic  loads  can  be  complex  to  define 
over  large  spatial  areas  and  a  simplifying  assumption  is  therefore  made  in  current 
practice.  The  hydraulic  loading  conditions  are  assumed  to  be  fully  dependent  in  terms 
of  recurrence  interval  (return  period)  within  a  flood  area  (i.e.  the  multivariate  variable 
of  the  hydraulic  loading  is  reduced  to  a  univariate  distribution).  This  enables  the 
integration  of  the  density  of  the  hydraulic  loads  over  the  consequence  function  to  be 
undertaken  in  tenns  of  a  single  likelihood  of  hydraulic  load  using  a  simple  integration 
procedure: 


EADxY"  P 


rxj,+xi  X,  ^ 

-Jxl — t-  <  X  <  — — J— 
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g(xf) 
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where  g(x  )  is  the  expected  economic  damage  for  the  hydraulic  load  x,  and  q  is  the 

total  number  of  loading  levels  (return  periods)  used  in  the  analysis.  Flood  areas  are 
typically  less  than  10  km2  and  Eqn.  2  is  evaluated  independently  for  each  flood  area. 
To  assess  the  risk  at  spatial  scales  larger  than  a  single  flood  area,  the  results  for  each 
flood  area  are  aggregated.  It  is  of  note  however,  that  important  infonnation  relating  to 
the  nature  of  flood  events  can  be  lost  when  considering  only  the  EAD.  This  approach 
does  not,  for  example,  enable  assessment  of  the  probability  and  consequences 
associated  with  a  single  flood  event  that  occurs  at  spatial  scales  larger  than  a  flood 
area,  when  there  is  partial  dependence  between  the  hydraulic  loading  conditions. 
Analysis  of  single  event  consequences  at  large  spatial  scales  (ie  bigger  than  a  single 
flood  area)  can  be  of  interest  and  importance  for  the  insurance  industry,  where  single 
event  loss  damages  are  of  interest  and  emergency  planners,  for  example,  where  large- 
scale  evacuation  plans  are  developed.  Further  discussion  on  these  issues  is  provided 
by  Environment  Agency  (2011b) 

Within  the  analysis  described  here,  the  simplifying  assumption  of  full  dependence 
within  and  between  flood  areas  is  replaced  with  a  multivariate  extreme  value  model. 
This  model  is  described  in  more  detail  below.  The  use  of  the  multivariate  approach 
enables  information  on  the  probability  and  consequences  associated  with  single  events 
to  be  captured  more  accurately,  whilst  also  considering  the  likelihood  of  the  flood 
defences  failing. 


Description  of  the  modelling  method 


Overview 

The  objective  of  the  new  modelling  system  is  to  improve  the  handling  of  the  spatial 
dependence  of  flood  events  within  the  system  risk  analysis  method  and  hence  to  solve 
Eqn.  1  with  an  explicit  consideration  of  the  dependence,  in  the  extremes,  of  the 
hydraulic  loading  conditions  (the  vector  A)  at  different  spatial  locations.  A  flow 
diagram  showing  an  overview  of  the  new  modelling  system  is  shown  in  Figure  1. 
Figure  2  shows  a  conceptual  diagram  of  the  primary  spatial  components  within  the 
modelling  system.  Each  stage  of  the  new  modelling  system  is  described  in  more 
detail  below. 


Fit  multivariate  statistical  model  to  hydraulic 
loads  over  the  entire  model  domain. 


Figure  1 


Figure  2 


Simulated 
hydraulic  loading 
events  from  the 
multivariate 
statistical  model 


Statistical  model  for  hydraulic  loads 

There  are  a  wide  range  of  multivariate  extreme  value  models,  Coles  and  Tawn  (1991), 
Joe  et  al.  (1992),  Coles  and  Tawn  (1994)  Ledford  and  Tawn  (1996),  for  example.  All 
of  these  models  have  restrictive  assumptions  relating  to  the  nature  of  dependence. 
More  specifically,  in  the  joint  tail  region,  one  variable  is  assumed  to  be  independent 


of  or  asymptotically  dependent  on  the  other  variables.  The  Heffernan  and  Tawn 
(2004)  approach  makes  no  such  assumptions  and  is  therefore  the  model  of  choice  for 
this  study.  Further  discussion  on  the  motivation  for  use  of  the  Heffernan  and  Tawn 
(2004)  method  in  the  context  of  flood  risk  analysis  is  provided  in  Environment 
Agency  (2011b). 


Let  A)  =  (Aj,  .  ..,Xj)t  be  a  time  series  of  flood  related  hydraulic  loading  conditions  at  a 
collection  of  locations  within  the  area  of  interest.  Before  dependencies  between 
locations  are  considered,  the  extreme  loads  at  each  location  are  studied  marginally. 
For  this  the  standard  peaks-over-threshold  approach  of  Davison  and  Smith  (1990),  is 
used:  cluster  maxima  are  identified  from  the  time  series  and  the  excesses  above  a 
suitably  high  threshold  are  fitted  to  the  generalised  Pareto  distribution  (GPD).  This 
defines  a  probability  model  for  large  values  of  the  variable  A): 


P{Xi>x\Xi>ui}  = 


1  +  6 


A  1 


for  x  >  Ui, 


(3) 


Where  the  subscript  i  denotes  the  spatial  location  of  the  data,  [:>,  >  0,  e  91  are  the 
GPD  parameters  and  [v]+  =  max(y,  0).  The  threshold  u,  is  chosen  to  be  just  large 
enough  to  ensure  a  stable  estimate  for  the  shape  parameter  <f,-  for  all  larger  thresholds. 


In  order  to  separate  the  marginal  characteristics  from  the  dependence  analysis,  it  is 
usual  to  standardise  the  data  to  common  margins  using  the  probability  integral 
transform.  This  is  often  referred  to  as  the  copula  approach.  The  Heffernan  and  Tawn 
(2004)  model  uses  standard  Gumbel  marginal  distributions  which  are  obtained  by 

setting  Yt  =  -log(-log[/f  (Af)])  where  Fj  is  an  estimate  of  the  cumulative 
distribution  function  for  A)-  For  this  the  GPD  fit  above  the  threshold  is  combined  with 
the  empirical  distribution  Ft  of  the  A)  values  to  give  the  following  semi-parametric 
function  (first  used  by  Coles  and  Tawn  (1991)): 
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The  transfonned  multivariate  time  series  Yt  =  (7,  Yd)t  retains  the  dependence 
structure  of  the  original  data  but  satisfies  P\Y!  >  y }  =  exp(-  exp[-  v])  for  each  different 
location  The  primary  aspect  of  the  Heffernan  and  Tawn  (2004)  approach  is  to  model 
the  dependence  between  extreme  values  of  7,  and  typical  values  of  the  remaining 
variables.  The  analysis  is  repeated  for  each  site  i  so  that  extreme  values  of  all 
variables  are  considered. 


Let  Y-i  denote  the  vector  of  all  variables  Yj  excluding  7,.  The  Heffernan  and  Tawn 
(2004)  approach  is  typically  applied  using  the  multivariate  non-linear  regression 
model 

Y-i  |  7  =  o  7  +  Yi  Z  for  7  >  v,  (5) 

where  v  is  a  high  threshold  on  7,  a  e  [0,  1]  and  b  <  1  are  vectors  of  parameters  and  Z 
is  a  vector  of  residuals.  Vector  arithmetic  should  be  interpreted  component- wise  so 


that  each  Yj  is  modelled  as  a  function  of  7  using  parameters  a/|;-  and  bf  ,  and  residual 

ZJ\i- 


The  regression  parameters  and  bj,  are  estimated  using  maximum  likelihood  under 
the  temporary  assumption  that  Zj\,  follows  a  Nonnal  Distribution  with  unknown  mean 
and  variance.  This  fit  uses  all  pairs  (7,  YJ)  corresponding  to  cluster  maxima  of  T,  >  v 
to  be  consistent  with  the  marginal  GPD  fits  made  to  cluster  maxima  of  Xh  Heffeman 
and  Tawn  (2004)  show  that  asymptotically  7  >  v  is  statistically  independent  of  the 
residual  Z;-|,.  The  threshold  v  is  therefore  chosen  to  be  just  large  enough  for  this 
condition  to  hold.  Once  all  parameter  estimates  have  been  found  a  non-parametric 
estimate  of  the  joint  distribution  of  Z  is  constructed  from  the  empirical  distribution  of 
the  sample  residuals. 

The  above  description  assumes  variables  7,  and  Yj  occur  concurrently  so  is  not 
appropriate  for  modelling  temporally  dependent  data  with  extreme  events  that  are 
lagged  between  gauges.  Keef  et  al.  (2009)  overcome  this  deficiency  by  fitting  the 
conditional  model  of  7,/+r  |  7y  for  a  selection  of  lags  r  for  each  gauge  j  ±  i-  This 
allows  the  model  to  be  used  to  simulate  new  events  over  a  range  of  lags  so  that  the 
largest  values  in  each  time  window  need  not  occur  concurrently  and  this  approach  has 
been  adopted  here. 

Statistical  simulation  of  extreme  hydraulic  loading  events 

The  fitted  model  provides  parameter  estimates  aJiLT  and  r  for  all  locations  i,j  ±  i  and 
for  a  range  of  lags  r.  Additionally,  for  each  location  i,  an  empirical  sample  of  joint 
residuals  Z  is  available,  each  of  which  provides  values  zi\ «,  x  for  each  j  ±  i  and  lag  r. 
These  can  be  used  to  simulate  a  large  number  of  dependent  peak  events  Y,  each  of 
which  consists  of  a  single  peak  value  for  every  location  in  the  network  with  associated 
lags  between  peaks.  These  are  then  transformed  to  give  samples  of  spatially 
dependent  peak  flow  events  X  on  the  original  scale. 

To  simulate  an  event,  a  conditioning  site,  i  must  be  selected.  The  value  7,  is  sampled 
above  the  threshold  v  and  the  remaining  7,  values  are  sampled  from  the  fitted  model 
for  7_,  t  |  7.  In  order  to  control  the  proportion  of  events  where  each  site  is  most 
extreme,  the  value  7  is  constrained  to  be  largest  by  rejection  sampling.  Further 
discussion  on  this  approach  is  provided  in  Environment  Agency  (2011b). 

The  simulation  consists  of  repeating  the  following  steps,  after  selecting  a  conditioning 
site  i: 

1 .  Sample  a  value  7,  from  the  standard  Gumbel  distribution  conditioned  to 
exceed  v. 

2.  Independently  select  one  of  the  joint  residuals  Z  for  site  i. 

3.  Calculate  7  r  =  afii  JYi  +  Z  /|;  r  for  all  j  ±  i  and  lags  r. 

4.  For  each  j,  set  Yj  to  be  the  maximum  7;-,r  over  selected  lags  r. 

5.  The  joint  sample  Y  is  rejected  unless  7  is  maximum  over  all  gauges. 

These  simulated  events  have  Gumbel  marginal  distributions  which  can  be  transformed 
to  the  original  scales  using  the  probability  integral  transform,  inverting  the  transform 
applied  to  the  original  data  using  GPD  fits  and  the  empirical  distribution. 


To  select  the  proportion  of  events  where  each  gauge  is  most  extreme,  Lamb  et  al. 
(2010)  propose  simulation  from  the  fitted  model  without  rejection  and  this  approach 
has  been  adopted  here.  For  each  gauge  i,  a  large  sample  of  samples  is  used  to  estimate 
P{Yt  >  Yj\/j  ±  i  |  Yi  >  v},  the  probability  that  Y,  is  largest  when  i  is  the  conditioning 
gauge.  Since  all  simulated  events  have  at  least  one  threshold  exceedance,  this  can  be 
used  to  estimate  the  proportion  of  all  events  where  Y,  is  largest  and  hence  variable  X, 
is  most  extreme. 

The  output  from  this  analysis  is  a  simulated  set  of  hydraulic  loading  events.  Each 
event  comprises  a  peak  flood  variable  at  each  location.  An  example  of  fluvial  flow 
events  simulated  from  the  statistical  model  and  compared  to  the  observed  data  is 
given  on  Gumbel  scales  in  Figure  3.  The  yellow  points  are  cluster  maxima  in  Y,  from 
the  raw  time  series  above  the  threshold  v  (147  cluster  maxima  in  Yt  were  found  above 
v  =  4.15  (on  Gumbel  scale)  from  a  27  year  dataset).  These  were  fitted  (ignoring  lags) 
to  give  the  blue  samples,  constrained  to  7/  >  Yj. by  rejection.  Similarly,  the  green 
points  are  cluster  maxima  in  Yj  fitted  to  give  the  red  samples  for  Yj  >  Y,.  The  diagonal 
defines  the  threshold  for  the  rejection  during  the  simulation  (i.e.  for  the  blue  dots 
when  Yj  >  Y, ).  It  is  apparent  that  the  dependence  in  the  extremes  for  site  i  within  the 
simulated  data  reproduces  well  that  within  the  observations.  Further  examples  are 
available  within  Heffernan  and  Tawn  (2004)  and  for  flood  related  variables  within 
Lamb  et  al.  (2010)  and  Environment  Agency  (201  le),  for  example. 
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Figure  3 


Defence  system  state  reliability  and  floodplain  analysis 

The  approach  used  to  represent  the  perfonnance  of  the  flood  defences,  the 
propagation  of  floodwater  across  the  floodplain  and  the  evaluation  of  economic 
damages  is  described  by  Gouldby  et  al.  (2008a)  which  extends  the  earlier  work  of 
Hall  et  al.  (2003)  and  Environment  Agency  (2003).  A  summary  description  of  this 
method  and  how  it  has  been  integrated  with  the  multivariate  extreme  value  method  is 
provided  below. 

Each  hydraulic  loading  event  that  has  been  simulated  from  the  multivariate  statistical 
model  forms  the  boundary  conditions  for  the  flood  defence  system  reliability  analysis 
and  subsequent  inundation  and  economic  consequence  evaluation.  The  line  of  flood 
defences  (embankments,  for  example)  that  fonns  the  boundary  between  the  river  or 
coast  and  the  floodplain  area  is  discretised  into  sections  that  have  different 
perfonnance  characteristics  (i.e.  type  of  defence,  elevation  and  condition),  see  Figure 
2.  Fragility  curves  are  used  to  define  the  performance  of  each  defence  section, 

USACE  (1996),  Simin  et  al  (2009),  Schultz  et  al.  (2010),  Vorogushyn  et  al.  (2009), 
for  example,  and  the  perfonnance  of  each  defence  section  is  assumed  to  be 
independent  from  any  other.  It  is  assumed  that  each  defence  can  exist  in  two  possible 
states,  failed  (ie  structural  failure  or  breached),  or  undamaged.  During  a  flood  event, 
water  can  enter  the  floodplain  through  any  particular  defence  if  it  is  a  breached 
defence,  or  if  it  remains  undamaged  but  is  overtopped.  As  there  are  a  finite  number  of 
defence  sections  and  only  two  possible  system  states  for  each  section,  there  are  a 
finite  number  of  defence  system  states  (combinations  of  failed  and  undamaged 
defences).  For  each  hydraulic  loading  event  it  is,  in  principle,  necessary  to  simulate 
the  flood  wave  propagation  and  associated  economic  damages  for  each  (2")  system 
state.  It  is,  however,  computationally  impractical  to  run  an  inundation  simulation  for 
all  of  the  different  defence  system  states  and,  for  each  hydraulic  loading  event,  a 
Monte  Carlo  simulation  of  the  defence  system  states  is  therefore  implemented. 

For  each  simulated  defence  system  state,  floodplain  discharge  volumes  are  calculated 
for  each  flood  defence.  This  water  is  then  propagated  across  the  floodplain  using  a 
computationally  efficient  flood  propagation  algorithm,  the  Rapid  Flood  Spreading 
Model  (RFSM),  Gouldby  et  al.  (2008a),  Lhomine  et  al.  (2008).  The  algorithm  uses  a 
pre-process  that  analyses  the  natural  floodplain  topography  to  distinguish  a  series  of 
localised  depressions.  The  depressions,  termed  Impact  Zones,  are  irregular  in  shape 
and  comprise  any  number  of  regularly  spaced  Impact  Cells  (the  underlying 
topographical  elevation  model),  see  Figure  2.  The  geometrical  properties  of  the 
Impact  Zones  are  stored  in  database  tables  and  these  fonn  the  computational  elements 
for  the  inundation  model.  The  volume  of  water  discharged  over  (or  through  breached 
defences),  is  then  spread  across  the  floodplain,  using  the  stored  infonnation  on  the 
Impact  Zones  and  a  final  flood  extent  and  floodplain  depths  for  each  Impact  Cell 
evaluated. 

As  the  model  is  volume  based,  it  is  exceptionally  computationally  efficient,  at  the 
expense  of  complete  representation  of  the  underlying  physical  processes.  As  there  is 
no  temporal  evolution  of  the  flood  wave,  it  is  not  possible  to  compute  velocities,  for 
example.  Recent  research,  Jamieson  et  al.  (2012),  has  however,  seen  the  development 
of  a  new  dynamic  model,  RFSM  EDA,  that  utilises  the  same  meshing  system  as  the 
RFSM  but  implements  a  new  simplified  formulation  of  the  shallow-water  equations, 
Bates  et  al.  (2010).  The  implementation  of  this  type  of  dynamic  model  within  the 


system  described  here,  affords  a  further  opportunity  for  significant  improvement.  It  is 
likely  the  type  of  steady-state  modelling  system  described  here  will  evolve  into  a 
more  dynamic  system  in  the  future.  In  principle,  it  would  also  be  possible  to  utilise 
the  infonnation  on  temporal  lags,  generated  within  the  application  of  the  multivariate 
analysis,  to  generate  time  varying  hydraulic  boundary  conditions. 

Flood  depths  output  for  each  Impact  Cell  are  then  combined  with  information  from 
the  National  Property  Database  and  the  widely  applied  functions,  of  Penning-Rowsell 
et  al.  (2005)  to  determine  economic  damages. 

Within  any  particular  discrete  Impact  Cell  within  the  floodplain,  the  probability  of 
exceeding  any  particular  flood  depth  (h),  conditional  on  the  specified  hydraulic 
loading  level  (x),  is  approximated  by: 

P(H>h\x)*r/m)  (6) 

where  m  is  the  total  number  of  defence  system  state  Monte  Carlo  realisations  and  /«/, 
is  the  number  that  exceed  h,  under  the  specified  loading  level  x.  Similarly,  the 
conditional  probability  of  exceeding  a  specific  economic  damage  level  (c)  is  given  by: 

P(C>c\x)Mm/m)  (7) 

Combining  the  Monte  Carlo  realisations  for  the  hydraulic  loading  with  the  defence 

system  states,  it  then  follows  that  the  unconditional  annual  probability  of  exceeding 
any  particular  level  of  economic  damage  is: 

P(C  >c)~  (8) 

ny 

where  lc  is  the  number  of  Monte  Carlo  realisations  that  exceed  c,  /  is  the  total  number 
of  Monte  Carlo  realisations  and  ny  is  the  number  of  years  of  data  that  has  been 
simulated.  The  risk,  expressed  in  the  usual  terms  of  EAD  is: 

i>, 

EAD  «  7=i_  (9) 

ny 

where  c,  is  the  economic  damage  arising  on  the  /th  realisation  of  the  hydraulic  loading 
conditions.  The  number  of  realisations  required  can  be  controlled  in  a  number  of 
different  ways,  for  example,  through  specification  of  convergence  criteria  on  the 
quantity  of  interest. 


Case  study  application 
Study  area  description 

The  case  study  location  is  the  Eden  catchment  in  the  northwest  of  England  (Figure  4). 
The  Environment  Agency  (2009a)  provides  an  in  depth  discussion  of  the  flooding 
problems  and  ongoing  management  of  flooding  in  the  catchment  and  this  is  briefly 
summarised  here.  The  Eden  catchment  covers  approximately  2400km2,  with  the 
primary  watercouses  being  Eden,  Irthing,  Eamount,  Petteril  and  the  Caldew.  Over 
90%  of  the  catchment  is  rural  and  the  total  population  is  around  244,000.  The  upper 
catchment  (in  the  South),  is  dominated  by  high  ground  (fells).  In  the  Lower  Eden 
Valley,  the  floodplain  widens,  affording  storage  capacity  during  times  of  elevated 
water  levels.  In  the  higher  regions  to  the  South,  the  average  annual  rainfall  exceeds 
2800mm  and  is  over  3  times  the  national  average  of  England  and  Wales  (920mm). 


The  main  population  centres  are  Carlisle,  Penrith  and  Appleby.  The  floodplain 
protection  varies  significantly  across  the  study  area  and  includes  high  ground  in  rural 
areas,  embankments  and  a  range  of  heavily  engineered  vertical  structures  within  urban 
areas,  Carlisle  for  example.  The  study  area  has  been  separated  into  10016  flood  areas 
and  is  consistent  with  those  used  in  the  National  flood  risk  assessment  of  England, 
Environment  Agency  (2009b).  It  is  however,  of  note  that  the  underlying  data  has 
been  modified  to  a  certain  degree  and  the  results  shown  here  do  not  reflect  the 
absolute  magnitude  of  risk  in  the  catchment. 


Figure  4 


Carlisle  has  a  history  of  flooding  with  significant  events  having  been  recorded  in 
1963,  1968,  1979,  1980,  1984  and  most  recently  in  January  2005,  when  nearly  3000 
homes  were  flooded,  3  people  died  and  the  resulting  flood  losses  were  estimated  as 
£400m,  Geographical  Association  (2009).  The  flooding  arose  as  a  result  of  a 
prolonged  period  (1  month)  of  high  rainfall  in  the  vicinity  of  Carlisle,  saturating  the 
soil  and  raising  local  ground  water  levels.  This  was  then  followed  by  a  period  of 
intense  rainfall  (2  months-worth  in  2  days),  over  the  period  6th  -8th  of  January, 
Geographical  Association  (2009).  Figure  5  shows  the  hydrograph  on  the  River  Eden 
for  this  period. 
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Figure  5 


Hydraulic  load  analysis 

Synthetic  flow  data  form  the  basis  of  the  analysis  and  these  stem  from  a  series  of 
linked  models.  Time  series  rainfall,  on  a  25km  spatial  grid  and  an  hourly  time 
interval  provide  the  driving  data  for  a  gridded  hydrological  model  (Grid-to-Grid,  or 
“G2G”,  Bell  (2007)).  These  rainfall  data,  together  with  estimates  of  potential 
evaporation  (PE)  have  been  generated  from  the  UK  Met.  Office  Hadley  Centre  global 
climate  model  (GCM),  HadCM3,  dynamically  downscaled  to  25km  using  the 
HadRM3  regional  model,  Jones  R  et  al.  (  2004).  The  G2G  Model  operates  on  a  1  km 
grid  covering  the  whole  of  the  UK,  and,  using  a  gridded  representation  of  the 
kinematic  wave  equation,  simulates  surface  and  sub-surface  flows.  A  full  description 
of  the  G2G  model  is  provided  in  Bell  (2007). 

The  G2G  Model  provides  27  years  of  time  series  flows  on  a  1km  UK  grid  along  the 
river  network.  The  gridded  data  were  first  interpolated  to  each  of  41 19  nodes  that 
span  the  river  network  by  selecting  the  largest  value  of  each  surrounding  grid  point  at 
every  time  step.  To  reduce  the  computational  burden,  rather  than  applying  the 
statistical  analysis  to  every  node,  node  points  that  had  a  value  of  Pearson’s  Product 
Moment  Correlation  over  0.99  with  a  remaining  node  were  removed  to  leave  a 
representative  subset  of  248  sites.  This  threshold  correlation  value  was  considered 
sufficiently  accurate  for  use  on  this  case  study  but  it  is  possible  a  further  reduction  in 
the  number  of  sites  could  have  been  achieved  using  a  lower  value  without  loss  of 
accuracy. 


The  marginal  extremes  at  each  of  these  selected  nodes  have  been  analysed  separately. 
Peak  flow  clusters  were  identified  using  the  runs  method  of  Smith  and  Weissman 
(1994)  before  the  cluster  maxima  above  a  high  threshold  were  fitted  to  GPD.  After 
transformation  to  Gumbel  marginal  distributions,  cross-correlations  were  analysed 


using  node  pairs  lagged  by  up  to  ±6  days.  In  the  majority  of  cases,  the  largest 
correlation  was  obtained  by  a  lag  within  ±2  days.  Hourly  lags  of  up  to  ±3  days  were 
therefore  used  to  fit  the  multivariate  statistical  model  to  every  pair  of  nodes  in  the 
subset,  adopting  a  similar  approach  to  that  of  Lamb  et  al.  (2010). 

Peak  flows  were  then  paired  with  local  maxima  at  every  other  node  so  that  extreme 
events  that  had  at  least  one  threshold  exceedance  could  be  identified.  709  such  events 
were  counted  giving  an  average  of  26  extreme  events  per  year.  With  this,  1000  years 
worth  of  hydraulic  loading  conditions  were  simulated  from  the  fitted  model  on  the 
Gumbel  scale.  Each  sample  comprised  a  peak  flow  level  for  all  248  selected  nodes 
and  had  at  least  one  threshold  exceedance.  Rejection  sampling  was  used  to  control  the 
proportion  of  events  that  are  most  extreme  for  each  site. 

After  transformation  back  to  the  original  scale,  linear  relationships  identified  between 
nodes  with  large  correlations  (>0.99)  were  used  to  extend  the  peak  flow  samples  to  all 
4119  river  nodes.  These  were  then  converted  to  water  levels  using  the  Conveyance 
Estimation  System,  Me  Gahey  et  al.  (2008)  and  mapped  to  all  34006  discrete  river 
defence  sections  in  the  catchment. 

Figure  6  shows  a  comparison  of  hydraulic  loading  data  (peak  fluvial  flows)  simulated 
from  the  statistical  model  for  a  subset  of  stations.  The  figure  shows  the  range  of 
dependence  in  the  hydraulic  loading  conditions  that  is  present  across  the  study  area 
and  highlights  the  flexibility  (in  terms  of  the  range  of  dependence  in  the  extremes) 
required  of  the  multivariate  model  to  capture  this.  Perhaps  unsurprisingly,  sites  in 
close  proximity  exhibit  a  high  degree  of  dependence  whilst  others,  further  a-field  are 
only  partially  dependent  on  one  another. 


Figure  6 
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Figure  7  shows  the  variation  in  return  period  of  the  hydraulic  loading  levels  across  the 
study  area  for  a  subset  of  single  realisations  from  the  multivariate  model.  The  spatial 
variation  in  the  intensity  of  the  flood  events  is  apparent.  Some  events  are  spatially 
diverse  and  extreme  flows  arise  across  the  entire  study  area,  whereas  other  events  are 
very  much  localised  with  extremes  arising  in  relatively  small  areas  over  the 
catchment.  The  event  shown  in  the  lower  left  of  Figure  7  is  of  particular  interest.  The 
event  comprises  relatively  high  flows  in  the  upstream  areas  of  the  catchment 
(generally  >  200  year  return  period),  to  the  South  but  the  downstream  areas,  in  the 
vicinity  of  Carlisle,  to  the  North,  flows  are  significantly  less  (generally  <50  year 
return  period).  This  highlights  the  complexity  of  the  hydrological  regime  in  the 
catchment  and  gives  an  indication  of  the  high  influence  of  localised  rainfall,  in  the 
lower  catchment,  as  evidenced  by  the  2005  extreme  event  (Figure  5). 


The  range  of  flood  events  that  are  shown  highlights  the  communication  difficulties 
that  have  arisen  in  the  past  when  describing  flood  events  in  terms  of  the  severity  of 
their  loading  conditions  -  “The  event  was  a  1  in  100  year  (flow)  event”,  for  example. 

It  is  apparent  that  statements  about  flood  events,  expressed  in  terms  of  the  extremes  of 
the  hydraulic  loading  conditions,  are  often  only  valid  for  small  spatial  areas. 
Attempting  to  express  widespread  flood  events,  like  those  that  affected  England  and 
Wales  in  the  summer  of  2007,  in  terms  of  return  period  of  the  hydraulic  loading 


conditions,  as  often  seems  to  be  a  requirement,  can  result  in  misleading  information 
and  a  general  misunderstanding  of  the  nature  of  the  flood  event. 
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Figure  7 

Lewis  et  al.  (2011)  draw  attention  to  the  validity  of  making  an  assumption  of 
complete  dependence  of  hydraulic  loading  level  (equal  return  period)  within  a  flood 
area  with  system  based  flood  risk  analysis  models,  for  example.  This  assumption  is 
made  to  enable  a  simplified  integration  of  the  joint  density  of  the  hydraulic  loads  over 
the  consequence  function  (Eqn.  2).  Their  analysis  was  conducted  within  the  context 
of  coastal  systems  and  it  is  possible  to  explore  the  validity  of  this  assumption  further 
here,  within  the  context  of  this  fluvial  system.  Figure  8  shows  the  variation  of  return 
period  within  one  of  the  largest  flood  areas  for  the  same  events  shown  in  Figure  6. 

The  flood  area  comprises  30  defences,  the  maximum  and  minimum  return  periods  and 
average  return  period  for  each  of  the  hydraulic  loading  events  are  shown  in  Table  1. 
Based  upon  the  synthetic  data  analysed,  it  is  apparent  that  even  within  a  flood  area 
there  can  be  significant  variation  in  the  return  period  of  the  water  level  and  the 
assumption  of  equal  likelihood  can  potentially  introduce  significant  bias  in  some 
cases.  This  reinforces  the  observations  on  coastal  systems  of  Lewis  et  al.  (2011). 


Figure  8 


Return 

period 

flow 

(years) 


.05 

.1 

.2 

.5 

1 

2 

5 

10 

20 

50 

100 

200 

500 

1000 

2000 


Flood 

area 


Economic  damages  and  risk  estimation 

Within  the  current  National  Flood  Risk  Assessment  of  England  and  Wales, 
Environment  Agency  (2009b),  the  EAD  is  derived  at  the  level  of  a  flood  area  and  then 
aggregated  to  obtain  a  regional  or  national  value  of  risk.  The  method  does  not  enable 
estimates  of  return  period  damages,  for  single  events  that  are  greater  than  a  flood  area. 
The  introduction  of  the  multivariate  method  does  however,  enable  this.  For  each 
realisation  of  the  multivariate  model  of  the  hydraulic  loads,  a  Monte  Carlo  simulation 
of  the  defence  system  states  and  associated  inundation  and  consequence  analysis, 
within  each  flood  area,  has  been  undertaken.  The  resulting  risk,  in  terms  of  EAD, 
over  the  whole  study  area  is  £19.7  million.  Figure  9  shows  the  spatial  variation  in  the 
risk.  Whilst  damages  occur  throughout  the  catchment,  the  greatest  concentration  of 
risk  is  in  the  north  in  the  vicinity  of  Carisle,  the  site  of  recent  and  historic  flooding. 


Figure  9 
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The  distribution  of  damage  aggregated  over  the  whole  study  area  is  shown  in  Figure 
10.  A  comparison  has  also  been  made  under  the  assumption  of  full  dependence  and 
complete  independence  of  hydraulic  load,  on  a  flood  area  basis.  It  is  interesting  to 
note  that  in  this  particular  study  area  the  results  show  the  modelled  dependence  of 
economic  damages  to  be  close  to  full  dependence.  This  is  in  contrast  to  the  extreme 
hydraulic  loading  conditions  that  showed  a  wide  range  of  dependence  across  the  study 
area.  The  high  dependence  in  damages  results  from  the  localised  nature  of  some 
extreme  rainfall  events  in  the  North  of  the  catchment,  causing  major  damage,  as 
evidenced  by  the  event  in  January  2005.  It  is  important  to  note  that  this  finding  is 
only  applicable  for  this  study  area  and  other  study  areas  will  differ.  For  example, 
Lamb  et  al.  (2010)  found  the  degree  of  dependence  was  closer  to  independence  for  a 
study  area  in  the  northeast  of  England. 


Validation 

Validation  is  an  important  consideration  for  any  numerical  model.  Risk  however,  is 
an  abstract  quantity  that  cannot  be  measured  and  hence  one  of  the  usual  routes  of 
verifying  numerical  models  through  comparison  with  measured  data  is  not  possible. 
Some  confidence  in  the  new  approach  described  here  can  be  gained  from  previous 
analysis.  For  example,  the  systems  risk  model  has  been  verified  for  use  by  the 
Environment  Agency  on  the  major  Thames  Estuary  2100  study,  Gouldby  et  al. 

(2008b)  and  the  National  flood  risk  assessment  of  England  and  Wales  Environment 
Agency  (2009b).  Lamb  et  al.  (2010)  and  Environment  Agency  (2011b),  provide  an  in 
depth  exploration  of  the  application  of  the  multivariate  method  of  Heffernan  and 
Tawn  (2004)  within  the  context  of  flooding.  Further  work  is  however,  required  with 
respect  to  the  validation  of  probabilistic  flood  risk  models.  Environment  Agency 
(2010b)  describes  a  framework  for  undertaking  this  type  of  analysis.  In  particular,  the 
framework  recognises  the  need  to  undertake  uncertainty  analysis  that  includes  both 
model  structural  as  well  as  input  data  and  model  parameter  uncertainties.  Whilst  prior 
work  on  the  propagation  of  input  data  and  parameter  uncertainty  has  had  some 
success,  Merz  and  Thieken  (2009a),  Gouldby  et  al.  (2010),  further  work  on  the 
handling  of  model  structural  uncertainties  is  required. 

There  are  a  number  of  possibilities  with  regard  to  assessing  model  structural 
uncertainties.  One  promising  approach  relates  to  the  development  of  benchmark  data 
sets  that  can  be  used  to  explore  and  verify  the  model  components  or  the  whole 
modelling  system.  For  example,  Environment  Agency  (2010a)  describes  a  range  of 
benchmark  tests  for  2D  inundation  models,  these  type  of  tests  can  be  used  to  gain 
insight  on  the  performance  of  reduced  complexity  inundation  models  (eg.  RFSM)  that 
are  used  within  the  probabilistic  system  models.  It  is  of  note  however,  this  type  of 
component  testing  does  not  necessarily  capture  the  full  range  of  model  structural 
uncertainties.  Apel  et  al.  (2009),  for  example,  describe  the  influence  of  upstream 
dyke  breaches  on  downstream  water  levels  and  hence  downstream  breach  likelihood. 


To  assess  this  type  of  model  structural  uncertainty,  integrated  river  channel/floodplain 
models  that  incorporate  probabilistic  defence  perfonnance  and  dynamic  breach 
growth  are  required,  Lhomine  et  al.  (2010),  for  example.  The  development  of  these 
benchmarking  data  sets  is  the  subject  of  ongoing  research  which,  when  complete,  is 
likely  to  afford  more  opportunities  for  validating  the  type  of  modelling  system 
described  here. 


Conclusions 

Flooding  is  a  global  problem.  Effective  flood  risk  management  decision  making  is 
underpinned  by  quantitative  system  models  of  flood  risk.  The  models  that  are  applied 
in  current  practice  within  England  and  Wales  incorporate  a  wide  range  of  hydraulic 
loading  conditions,  a  representation  of  the  likelihood  of  failure  of  the  flood  defence 
assets  and  a  computationally  efficient  inundation  model.  The  current  models  do 
however,  utilise  a  simplifying  assumption  to  facilitate  the  evaluation  of  the  flood  risk. 
The  hydraulic  loading  events  within  a  flood  area  are  assumed  to  be  fully  dependent  in 
terms  of  likelihood.  This  assumption  is  made  to  facilitate  the  calculation  of  risk  in  a 
computationally  efficient  manner  and  to  overcome  the  complexities  relating  to  the 
dependence  structure  within  the  extremes  of  the  hydraulic  loading  conditions. 

An  artefact  of  this  assumption  is  the  inability  of  the  current  method  to  provide 
information  on  the  likelihood  of  damages  for  single  flood  events  at  spatial  scales 
larger  than  a  flood  area.  Within  the  current  method,  risk,  expressed  as  EAD  is 
calculated  on  a  flood  area  basis  and  aggregated  to  obtain  risk  at  larger  spatial  scales. 
Infonnation  on  the  probability  of  exceeding  a  specified  economic  damage  threshold 
for  a  single  event  is  not  currently  available.  This  infonnation  can  be  important  for  the 
insurance  industry  and  emergency  planners,  for  example.  Recent  research  identified 
the  potential  of  a  multivariate  extreme  value  statistical  model  to  improve  this  aspect 
of  the  existing  modelling  system. 

To  overcome  these  deficiencies  within  the  current  method,  a  new  modelling  system 
has  been  developed.  The  new  system  integrates  a  multivariate  extreme  value  method 
of  the  hydraulic  loading  conditions  with  the  existing  modles  of  flood  defence 
performance,  inundation  and  economic  damages.  The  new  modelling  system  has 
been  applied  on  a  case  study  site  to  investigate  fluvial  flooding  in  a  catchment  in  the 
North  West  of  England.  A  Monte  Carlo  simulation  of  the  fitted  multivariate  statistical 
model  enables  the  hydraulic  boundary  conditions  across  the  flood  defence  system  to 
be  represented  more  realistically.  The  resulting  method  retains  the  system 
characteristics  of  the  current  methodology  and  its  ability  to  reflect  the  performance  of 
the  flood  defence  infrastructure.  This  is  implemented  through  a  secondary  Monte 
Carlo  simulation  of  the  defence  system  states,  conditional  on  the  output  of  the 
multivariate  model. 

This  development  has  enabled  the  exploration  of  the  assumption  of  full  dependence  of 
hydraulic  loading  likelihood  within  a  flood  area.  The  analysis  shows  that  in  some 
cases,  for  a  single  flood  event  the  return  period  can  vary  significantly,  even  within  a 
flood  area.  This  is  likely  to  introduce  a  bias  in  estimates  of  flood  risk  that  utilise  this 
assumption.  This  concurs  with  the  findings  of  recent  research  relating  to  coastal 
hydraulic  loading  events. 


The  implementation  of  the  multivariate  method  enables  the  probability  of  single  flood 
event  damages  to  be  quantified.  This  refinement  can  potentially  offer  insights  into  the 
spatial  characteristics  of  single  flood  events  that  are  not  currently  possible  with  the 
existing  MDSF2  system.  It  also  highlights  the  limitations  of  using  the  EAD.  EAD  is  a 
relatively  simplistic  measure  of  flood  risk  and  analysis  of  additional  infonnation 
comprised  within  the  full  distribution  can  be  of  importance. 

Some  confidence  in  the  system  developed  here  can  be  afforded  through  previous 
application  and  verification  of  the  different  modelling  components.  Validation  of 
system-based  probabilistic  models  remains  a  challenge.  Ongoing  research  on  the 
development  of  benchmark  tests  for  probabilistic  models  is,  however,  likely  to  afford 
greater  opportunity  for  rigorous  validation  in  the  future. 
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Figure  captions 


Figure  1  Flow  diagram  of  the  new  multivariate  system  based  modelling  system. 

Figure  2  Spatial  components  of  the  new  multivariate  system  based  model. 

Figure  3  Comparison  of  observed  peak  flow  events  (cluster  maxima  in  Yt  and  Yj, 
yellow  and  green  respectively)  against  samples  from  the  fitted  statistical  model 
(ignoring  time  lags)  for  conditioning  site  i  and  j  (blue  and  red  respectively)  plotted  on 
the  Gumbel  scale. 

Figure  4  Map  showing  the  study  area  for  the  case  study  with  the  representative  nodes 
spanning  the  river  network.  Contains  Ordnance  Survey  data  ©  Crown  copyright  and 
database  right  2012. 

Figure  5  Hydrograph  for  the  River  Eden  showing  elevated  flows  several  weeks  prior 
to  the  January  2005  event  (data  source.  National  River  Flow  Archive). 

Figure  6  Comparison  of  peak  flow  events  sampled  from  the  statistical  model  for  four 
sites  on  the  river  network. 

Figure  7  Variation  in  return  period  water  levels  across  the  study  area  for  a  subset  of 
sampled  hydraulic  loading  events. 

Figure  8  Variation  in  return  period  water  levels  within  a  particular  flood  area. 

Figure  9  Spatial  variation  in  EAD  over  the  study  area. 

Figure  10  Return  period  total  damage  for  the  study  area  for  the  fully  dependent 
(green),  partially  dependent  using  the  statistical  model  (blue)  and  independent  (red) 
cases. 


Tables 


Table  1  Minimum,  mean  and  maximum  return  period  flows  over  the  30  defences  in  the  highlighted 
flood  area  for  each  of  the  events  in  Figure  5.  The  events  in  the  figure  are  numbered  in  rows  beginning 
in  the  top-left. 
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The  need  for  large-  (region-)scale  probabilistic  simulations  means  that  2D  inundation  models  are  still  limited  by 
computational  requirements.  In  addition  to  parallelisation  and  physical  process  simplification,  attempts  to  reduce 
runtimes  typically  involve  coarsening  the  computational  mesh,  which  can  smooth  important  topographic  features 
and  hence  limit  accuracy.  This  paper  presents  a  new  2D  flow  model  that  uses  an  enhanced  diffusion-wave,  and 
incorporates  sub-element  topography  in  a  computational  mesh  that  adapts  to  the  terrain  features.  The  model  utilises 
a  fine  topographic  resolution  without  having  to  use  a  fine  computation  mesh,  and  so  achieves  fast  computational 
runtimes.  The  model  has  been  tested  against  the  Environment  Agency's  2D  benchmarking  tests,  and  even  though  the 
model  is  designed  to  operate  at  larger  spatial  scales  than  those  in  the  benchmarking  tests,  it  is  shown  to  provide 
comparable  accuracy  relative  to  a  selection  of  conventional  2D  models,  at  significantly  faster  computational  speeds. 
The  model  therefore  has  the  potential  to  offer  a  step  change  in  performance  of  large-scale  probabilistic  flood 
mapping  and  systems  flood  risk  analysis  modelling. 


Notation 

A i  impact  zone  plan  area  (m2) 

Ap  panel  flow  area  (m2) 

c  celerity  of  a  wave  (m/s) 

/  an  interface  between  impact  zones 

g  gravitational  acceleration  (m/s2) 

hj  impact  zone  water  depth  (m) 

i,  j  impact  zones 

n  Manning’s  coefficient  of  friction 

itj  directional  unit  vector  between  an  impact  zone  and  its 
neighbour  j 

Pp  panel  wetted  perimeter  (m) 

p  interface  panels/sections 

Qf  interface  flow  rate  (m3/s) 

Qp  panel  flow  rate  (m3/s) 

Rp  panel  hydraulic  radius  (m) 

Sf  interface  water  surface  slope  (m/m) 

t  time  (s) 

Ui  impact  zone  velocity  vector  (m/s) 

Vi  impact  zone  volume  (m3) 

w f  interface  width  (m) 

Xj  impact  zone  cross-sectional  flow  area  (m2) 

>/ f  interface  water  level  (m) 


rji  impact  zone  water  level  (m) 
a  constant  used  for  scaling  the  Courant  number 

[A  constant  used  for  velocity  calculation 

gp  frictional  wetted  height  on  panel  sides  (m) 

Ax  sub-element  (or  panel)  cell  width  (m) 

1.  Introduction 

Large-scale  flood  mapping  is  a  primary  requirement  of  the  Floods 
Directive  (EC,  2007),  and  probabilistic  flood  risk  models  that 
require  computationally  efficient  2D  components  are  in  increasing 
demand.  The  Environment  Agency  of  England  and  Wales’  (EA) 
national  flood  risk  assessment  (NaFRA)  (Environment  Agency, 
2009)  and  modelling  decision  support  framework  (MDSF2) 
(Environment  Agency,  2011)  have  utilised  simplified  inundation 
models  for  almost  a  decade,  due  to  the  number  of  simulations 
required  to  undertake  comprehensive  risk  analyses  (Gouldby  et 
al.,  2008a;  Hall  et  al.,  2003).  These  risk  models  have  also  been 
successfully  applied  for  a  wide  range  of  other  purposes  (e.g. 
Evans  et  al.,  2006;  Gouldby  et  al.,  2008b;  jWoodward  et  al.. 
Dll),  and  there  is  increasing  demand  to  improve  the  reliability 
of  the  results,  particularly  the  inundation  aspects  (National  Audit 
Office,  2011). 
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It  is  well  established  that  for  a  given  hydrological  input,  ground 
elevation  and  topographic  features  dominate  the  hydraulic  in¬ 
undation  process  (Romanowicz  and  Beven,  2003;  Zhang  and 
Cundy,  1989).  In  small-scale  studies,  particularly  urban  environ¬ 
ments,  computational  grid  sizes  must  be  of  the  order  of  1-5  m  to 
appropriately  characterise  the  underlying  topography  (Mark  et  al. , 
2004).  In  practice,  however,  large-scale  and  probabilistic  simula¬ 
tions  can  rarely  be  completed  at  this  resolution.  Reducing  the 
grid  size,  for  example,  has  a  dramatic  effect  on  the  computational 
cost  associated  with  full  shallow-water  equation  (SWE)  models 
and,  counterintuitively,  regular  grid  diffusion-wave  models  tend  to 
be  even  slower  at  such  resolutions  (Hunter  et  al,  2008). 

While  it  is  evident  that  2D  inundation  simulations  over  large 
areas  cannot  be  achieved  with  grids  of  equivalent  length-scale  to 
natural  topographic  variation,  using  traditional  grids  with  coarse 
resolution  can  artificially  smooth  important  topographic  features. 
To  address  this  shortcoming,  there  has  been  increasing  develop¬ 
ment  of  models  that  employ  a  sub-grid  representation  with  the 
aim  of  improving  topographic  detail  while  maintaining  computa¬ 
tional  efficiency  (Casulli  and  Stelling,  2011;  Hartnack  et  al., 
2009;  McMillan  and  Brasington,  2007;  Yu  and  Lane,  2006b).  Yu 
and  Lane  (2011)  found  that  post-processing  of  the  DEM  to  re¬ 
introduce  the  topographic  features  improves  the  simulation 
accuracy  of  their  sub-grid  approach.  This  demonstrates  that 
although  accounting  for  sub-grid  mass  storage  effects,  most  sub¬ 
grid  approaches  have  difficulty  representing  sub-grid  flow  block¬ 
age  effects  unless  the  grid  cell  boundaries  are  perfectly  aligned 
with  topographic  features.  An  exception  to  this  is  the  multi¬ 
layered  approach  of  Chen  et  al.  (2008),  which  can  cope  with 
urban  flow  blockages  within  a  grid  cell.  Yet  this  is  limited  to 
simplistic  building  layouts  and  it  would  be  difficult  to  apply  to 
real  catchment  topography.  For  coarse  modelling  this  issue  was 
avoided  in  the  past  by  using  manual  delineation  of  flood  cells 
along  floodplain  features,  such  as  railway  embankments  and 
dykes  (Estrela  and  Quintas,  1994;  Romanowicz  et  al.,  1996; 
Zanobetti  et  al.,  1970),  but  manually  creating  such  grids  is 
subjective  and  expensive  and  cannot  be  practically  undertaken  at 
a  large  scale. 

This  grid  issue  has,  to  some  extent,  been  resolved  by  the  rapid 
flood  spreading  method  (RFSM)  (Gouldby  et  al.,  2008a;  HR 
Wallingford,  2006),  wherein  computational  elements,  known  as 
impact  zones,  are  automatically  defined  to  precisely  follow 
topographic  features.  In  addition,  the  sub-element  topography  is 
resolved  at  the  level  of  the  underlying  digital  elevation  model, 
and  so  always  captures  the  critical  topographic  crests  causing 
flow  blockage,  without  the  need  to  manually  re-introduce  these 
features. 

The  RFSM  was  first  implemented  with  a  simple  spreading 
algorithm  that  conserved  volume  but  did  not  represent  the 
temporal  evolution  of  the  flood  wave  (direct  RFSM).  This  was 
later  improved  with  a  version  that  attempted  to  account  for  the 
frictional  and  dynamic  effects  of  floodplain  propagation,  using  a 


simplified  approach  (Lhomme  et  al,  2009).  More  recently, 
attempts  have  been  made  to  improve  process  representation  by 
incorporating  a  time-stepping  analytical  approximation  to  the 
diffusion  wave  (dynamic  RFSM)  that  is  similar  in  dynamics  to 
the  raster  storage  model  Lisflood-FP  (Bates  and  De  Roo,  2000). 
Output  from  this  model  compared  well  with  that  from  a  full  SWE 
model  on  a  large-scale  site  in  Ireland  (Lhomme  et  al.,  2012),  but 
less  well  when  employed  in  the  EA’s  2D  hydraulic  modelling 
benchmark  tests  (Wright  et  al,  2012).  Under  these  tests,  accuracy 
was  constrained  due  to  the  use  of  a  constant  time  step  and  flow 
limiters,  as  has  been  demonstrated  in  other  diffusive-type  models 
(Hunter  et  al.,  2005). 

This  paper  presents  a  new  version  of  the  RFSM  model  that 
overcomes  some  of  the  limitations  noted  above.  This  new  model, 
RFSM-EDA  (RFSM  -  explicit  diffusion  wave  with  acceleration 
term)  -  follows  the  sub-element  impact  zone  approach  but  uses  a 
new  formulation,  similar  to  the  diffusion  wave  but  incorporating 
the  local  acceleration  term  of  the  Saint  Venant  equations  (Bates 
et  al,  2010).  An  adaptive  time  step  has  been  implemented,  and 
all  flow  limiters  have  been  removed.  The  effectiveness  of  the 
model  is  demonstrated  using  a  selection  of  the  EA’s  2D  hydraulic 
benchmark  tests  (Wright  et  al.,  2012). 

2.  Model 

RFSM-EDA  is  based  on  the  same  mesh  concept  as  the  direct 
RFSM  (Gouldby  et  al.,  2008a;  Lhomme  et  al,  2009)  and  the 
dynamic  RFSM  (Environment  Agency,  2010,  Lhomme  et  al., 
2012).  See  Figure  1  for  a  mesh  schematic.  It  incorporates  the 
following  primary  assumptions. 

■  The  domain  can  be  divided  up  into  discrete  and  hydraulically 
consistent  topographic  depressions,  called  impact  zones  (IZs). 

■  The  water  surface  elevation  within  each  IZ  is  constant. 

■  The  relationship  between  water  surface  elevation  and  volume 
in  an  IZ  can  be  defined  by  a  non-hysteretic  relationship. 

■  The  flow  rates  between  neighbouring  IZs  are  calculated 
linearly  across  the  interface  between  them,  independently  of 
other  neighbours. 

■  The  interface  can  be  characterised  by  a  level-width 
relationship,  where  the  width  is  assumed  to  increase  with 
increasing  level. 

2.1  Pre-processing  algorithm 

Before  computation  can  begin,  the  IZs  are  defined  through  a  pre¬ 
processing  algorithm.  In  a  first  stage,  IZs  are  delineated  around 
collections  of  cells  which,  following  the  line  of  greatest  slope, 
would  drain  to  the  same  topographic  low  point.  This  produces  IZ 
boundaries  defined  along  topographic  crests  and  high  points.  In  a 
second  stage,  the  original  IZs  are  modified  to  ensure  that  they  are 
above  a  certain  minimum  area,  and  that  the  interfaces  between 
them  are  above  a  minimum  communication  depth.  These  mod¬ 
ifications  are  controlled  by  user-defined  parameters;  appropriate 
values  will  vary  depending  on  the  landscape  and  the  DEM 
resolution.  Finally,  the  level-volume  and  neighbour  level- width 


2 


Water  Management 

Volume  165  Issue  WM1 


A  highly  efficient  2D  flood  model  with 
sub-element  topography 

Jamieson,  Lhomme,  Wright  and  Gouldby 


Figure  1.  Schematic  of  an  impact  zone  with  a  neighbour,  in  plan 
and  profile.  Showing  irregular  boundaries  and  selected  key 
variables.  Solid  grey  represents  a  volume  of  water 


relationships  are  calculated,  and  the  results  are  written  to  a 
database. 

2.2  Governing  equations 

The  derivation  of  the  flow  equations  stems  from  the  approach  of 
Bates  et  al.  (2010).  Starting  with  the  one-dimensional  Saint 
Venant  equations,  advection  is  assumed  negligible  and  the  equa¬ 
tions  are  discretised  semi-implicitly,  but  rearranged  into  an 
explicit  form.  The  hydraulic  radius  is  calculated  in  full  including 
the  friction  on  the  side  of  cells.  This  differs  from  other  models 
(e.g.  Bates  et  al. ,  2010;  McMillan  and  Brasington,  2007;  Yu  and 
Lane,  2006a)  because  the  assumption  that  friction  is  only 
encountered  on  the  cell  base  may  not  be  appropriate  in  highly 
variable  terrain.  A  single  flow  is  required  for  each  neighbour 
interface,  /,  but  to  avoid  sudden  changes  in  hydraulic  radius  in 
complex  topography,  the  fluxes  are  evaluated  as  the  sum  of 
individual  fluxes  across  a  number  of  interface  panels,  equal  to  the 
number  of  sub-element  cells,/),  in  the  interface: 

w+A,  V-  (%-gAtAft) 

1.  f  ^l+gAtn^l/A^R'/'3 


where  (Qj  is  the  interface  flow,  Qp  is  the  panel  flow,  t  is  time,  g  is 
gravitational  acceleration,  Ap  is  panel  area,  Rp  is  hydraulic  radius 
of  the  panel,  n  is  Manning’s  coefficient  and  Sf  is  the  water 


surface  slope  across  the  interface.  See  Figures  1  and  2  for 
schematic  diagrams  showing  the  relationship  between  the  vari¬ 
ables. 


The  panel  hydraulic  radius  is  calculated  by 


The  panel  area,  Ap,  is  the  difference  between  the  interface  water 
level,  r]f,  and  the  panel  ground  level,  zp ,  multiplied  by  the  sub¬ 
element  cell  width,  Ax: 


3.  Alp  =  A x{r/f  -  zp) 


The  panel  wetted  perimeter,  Pp,  is  the  summation  of  the  wetted 
base  (i.e.  the  width  of  the  sub-element  cell)  and  the  wetted  height 
to  one  or  both  of  the  adjacent  panels,  gp: 

4.  Ptp  =  gp  +  Ax 


To  evaluate  Ap,  the  interface  flow  level,  rtf,  is  needed,  and  a 
number  of  different  approaches  can  be  applied.  Using  the  mean 
of  the  levels  in  the  adjacent  IZs  is  problematic  because  negative 
depths  occur  when  the  downstream  level  is  below  the  interface 
crest.  The  dynamic  RFSM  (Lhomme  et  al.,  2012)  avoids  this  by 
switching  to  the  upstream  level  when  the  mean  level  would  create 
a  negative  depth.  However,  this  can  cause  sudden  and  cyclical 
jumps  in  the  interface  depth  if  the  downstream  level  fluctuates 


Figure  2.  Schematic  of  an  interface  between  two  neighbouring 
impact  zones.  Solid  grey  colouring  represents  water  part¬ 
submerging  the  interface,  and  the  demarked  rectangle  represents 
a  calculation  panel  corresponding  with  an  individual  sub-element 
cell 
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around  the  level  of  the  interface.  In  RFSM-EDA  a  smoother 
result  is  obtained  by  always  using  the  upstream  level: 

5.  77}  =  max^',  77J) 


The  interface  slope,  Sf,  is  calculated  by  dividing  the  difference  in 
neighbouring  IZ  water  levels  by  the  separation  distance  between 
their  centroids. 

The  solution  is  progressed  by  applying  the  conservation  of  mass: 

V'+At  =  V[  +  At  J2Qf+A‘ 

6.  j 


where  Vt  is  the  volume  in  IZ  i,  and  j  is  an  IZ  neighbour  of  i.  Vi  is 
a  function  of  rjj,  the  IZ  water  level,  and  this  relationship  is 
defined  in  advance  in  look-up  tables  created  during  the  pre¬ 
processing  stage.  Therefore,  the  IZ  volume  can  be  efficiently 
converted  into  a  water  level  for  use  in  the  flux  calculations. 

2.3  Numerical  stability 

The  scheme  is  subject  to  the  Courant-Freidrichs-Lewy  (CFL) 
condition,  which  is  satisfied  by  ensuring  that  the  domain  of 
dependence  of  the  interfaces  of  an  IZ  should  not  exceed  the  area 
of  the  IZ,  as  used  by  Guinot  and  Soares-Frazao  (2006).  This 
version  of  the  CFL  condition  is  more  appropriate  for  irregular¬ 
shaped  elements  than  that  used  by  Bates  et  at.  (2010),  because  it 
uses  areas  rather  than  lengths.  It  also  differs  by  including  velocity 
with  celerity.  The  maximum  permissible  time  step,  Afmax,  is 
given  by 


Afmax  =  a  min 


7. 


^w/max(||«'|| 

j 


where  a  is  a  constant  used  to  scale  the  predicted  time  step,  A,  is 
the  surface  area  of  i,  Wf  is  the  interface  width,  iq  is  the 
magnitude  of  the  IZ  velocity  vector,  and  the  celerity  of  a  wave,  c, 
is  given  by 


interface  level  of  its  neighbours.  Similarly,  as  an  IZ  dries, 
inappropriately  large  flows  will  not  cause  a  negative  depth,  as  the 
stored  volume  below  the  minimum  interface  level  will  absorb  the 
excess  flux.  These  effects,  resulting  from  the  IZ  shape,  provide  a 
natural  resistance  to  model  instability.  Therefore  no  special 
wetting  or  drying  treatments  are  explicitly  represented  within  the 
model. 

2.5  Velocities 

The  velocities  calculated  at  the  interfaces  could  be  used  as  a 
surrogate  for  the  IZ  average  velocity  in  flat  topography,  but  this  is 
not  appropriate  when  the  IZs  have  a  depression-like  shape.  In  this 
case  the  interface  velocities  are  expected  to  be  relatively  shallow 
and  fast,  compared  with  deeper  and  slower  flow  conditions  at  the 
IZ  centre.  To  convert  the  interface  velocities  to  an  area-average 
velocity  vector,  an  additional  step  is  necessary.  Assuming  the  IZs 
are  of  regular  shape,  the  volume  of  water  that  has  been  fluxed  out 
of  the  IZ  (using  the  results  of  Equation  1)  is  divided  by  the  area 
of  a  representative  cross-section  through  the  centre  of  the  IZ,  Xf. 


E  Q'jnj 

- J— -  where  Q'  >  0 


where  «7-  is  the  unit  vector  between  the  IZ  and  neighbour 
centroids,  used  to  provide  the  velocity  as  a  vector.  Whether  the 
IZ  shape  is  assumed  cubic,  cylindrical  or  as  an  inverted  cone,  the 
calculation  for  the  IZ  cross-sectional  area  can  be  written  as 

10.  x\  =  pJ(W) 


where  /3  is  a  constant  which  for  the  aforementioned  shapes  takes 
on  a  value  between  0-96  and  IT 3.  As  we  assume  the  IZs  to  be  of 
variable  shapes  and  sizes,  /3  is  given  a  value  of  1  for  simplicity.  It 
is  important  to  note  that  this  velocity  does  not  impact  on  the 
fluxes  between  IZs,  which  are  calculated  independently.  The  only 
impact  it  has  on  the  model  is  through  the  CFL  condition 
(Equation  7). 


8.  c'i  =  \[gh'i 


where  hi  is  the  depth  of  water  in  IZ  i. 

2.4  Wetting/drying 

In  some  reduced  complexity  models  an  algorithm  is  used  to 
reduce  over-rapid  wetting  or  drying  (Bradbrook  et  at.,  2004;  Yu 
and  Lane,  2006a).  As  the  IZs  are  assumed  to  have  topographic 
barriers  as  crests  between  them,  when  an  IZ  initially  wets,  the 
water  cannot  leave  until  it  fills  the  volume  below  the  lowest 


3.  Application 

3.1  ■■■■■■ 

3.1.1  Environment  Agency  benchmark  tests 

The  EA  has  produced  a  set  of  hydraulic  benchmark  tests  designed 
to  test  a  range  of  predictive  abilities  of  2D  inundation  models. 
Details  of  the  test  specifications  are  provided  in  Wright  et  at. 
(2012)  and  Environment  Agency  (2010),  so  they  are  only  briefly 
described  in  this  paper.  RFSM-EDA  has  been  assessed  on  most 
of  these  tests,  although  only  the  results  of  tests  2A,  4,  5  and  8A 
are  shown  here.  Table  1  provides  a  breakdown  of  the  tests,  with  a 
justification  provided  for  those  not  shown. 
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Test 

Shown  in  paper? 

Reason 

1 

X 

RFSM-EDA,  as  most  models,  performs  well  on  this  test,  but  there  is  little  of  interest  in  the  results 

2A 

/ 

See  Section  3-2 

3 

X 

A  test  of  momentum  conservation,  which  the  numerical  scheme  of  RFSM-EDA  is  not  expected  to 
achieve 

4 

/ 

See  Section  3-3 

5 

/ 

See  Section  3-4 

6 

X 

Dam  break  scenarios  require  a  full  SWE  scheme  with  shock-capturing  ability,  so  RFSM-EDA  is  not 
tested 

7 

X 

Dynamic  linking  with  a  ID  element  has  not  yet  been  tested  for  RFSM-EDA 

8A 

/ 

See  Section  3-5 

8B 

X 

Dynamic  linking  with  a  ID  element  has  not  yet  been  tested  for  RFSM-EDA 

Table  1.  Benchmarking  tests  that  have  been  completed,  and 
justification  for  those  not 

3.1.2  Comparison  with  other  models 

RFSM-EDA  has  been  compared  against  a  number  of  other 
models  to  provide  a  context  for  the  results,  rather  than  to  draw 
specific  conclusions  about  these  individual  models.  While  this  is 
not  a  rigorous  test  of  the  model’s  validity,  in  the  absence  of 
validation  data  the  model  is  compared  with  a  range  of  respected 
and  widely  used  models.  Two  finite-volume  SWE  models  are 
shown,  InfoWorks-ICM  (Innovyze,  2011;  Lhomme  et  al.,  2010) 
and  Tuflow-FV  (2nd-order  spatial  accuracy)  (Environment 
Agency,  2010).  Three  simplified  models  are  also  shown:  JFLOW- 
GPU,  a  regular  grid  diffusion- wave  model  (Bradbrook  et  al., 
2004;  Lamb  et  al.,  2009);  the  dynamic  RFSM,  also  a  diffusion- 
wave  model  but  with  the  same  sub-element  representation  as 
RFSM-EDA  (Environment  Agency,  2010;  Lhomme  et  al.,  2012); 
and  Lisflood-ACC,  which  has  a  similar  numerical  approach  to 
RFSM-EDA  but  is  based  on  a  regular  grid  (Bates  et  al.,  2010; 
Neal  et  al.,  2011).  For  these  tests  all  the  models  adhered  to  the 
test  specifications  apart  from  the  dynamic  RFSM,  which  used  an 
equivalent  (though  not  identical)  mesh  to  ‘mesh  A’  used  by 
RFSM-EDA,  described  in  the  following  section. 

3.1.3  Application  of  RFSM-EDA 

The  primary  results  for  RFSM-EDA  are  created  using  a  mesh 
significantly  coarser  than  in  the  other  models,  but  with  a  sub¬ 
element  cell  resolution  corresponding  to  the  specified  grid 
resolution  of  the  tests.  This  is  called  mesh  A.  However,  some 
extra  simulations  have  been  carried  out  using  different  computa¬ 
tional  meshes  that  offer  additional  insight.  Mesh  B  uses  a 
similarly  coarse  computational  grid,  but  utilises  the  finest  topo¬ 
graphic  resolution  available  in  the  raw  DTM  for  its  sub-element 
resolution.  Mesh  C  replicates  the  test  specification  exactly,  like 
the  other  models.  This  means  using  a  fine  computational  mesh, 
with  each  RFSM-EDA  mesh  element  containing  one  topographic 
sub-element  cell.  A  summary  of  the  three  mesh  types  for  the 
different  tests  is  provided  in  Table  2. 

For  meshes  A  and  B,  the  results  have  been  produced  with 


significantly  coarser  meshes  than  recommended,  and  this  should 
be  noted  when  considering  the  results.  For  example,  the  compara¬ 
tive  models  have  extracted  results  from  small  grid  cells  contain¬ 
ing  only  the  specified  test  points,  whereas  RFSM-EDA  uses 
considerably  larger  computational  elements,  which  may  represent 
the  hydraulic  conditions  not  just  in  the  location  of  the  test  points 
but  at  distal  locations  as  well. 

Mesh  C  has  been  used  for  tests  2A  and  5  for  comparative 
purposes,  but  in  practice  RFSM-EDA  would  not  be  used  on  such 
a  mesh,  as  there  are  no  benefits  in  using  the  IZ  methodology 
when  each  IZ  contains  only  one  sub-element  cell.  In  fact,  the 
additional  computational  overhead  of  the  sub-element  approach 
(e.g.  calling  volume/level  look-up  tables)  makes  the  use  of  IZs 
with  one  sub-element  cell,  or  only  a  few  cells,  more  costly  than 
using  ‘traditional’  grids. 

All  the  RFSM-EDA  simulations  were  completed  on  a  machine 
running  Windows  XP  with  a  3  0  GHz  processor  and  8  GB  of 
RAM,  connecting  to  an  SQL  database  on  a  network  server. 

3.2  Test  2:  Filling  of  floodplain  depressions 
Test  2  is  designed  to  demonstrate  a  model’s  ability  to  deal  with 
inundation  processes  in  a  low-momentum  event.  The  test  is  a 
square  domain  of  16  topographic  depressions,  with  test  points  in 
each,  and  an  inflow  hydrograph  in  the  top  left  corner.  This  is  an 
extreme  test  of  the  IZ  schematisation;  rather  than  the  specified 
—  10  000  elements,  the  RFSM-EDA  mesh  A  uses  only  16 
automatically  generated  elements,  one  per  depression.  This  means 
that  mesh  A  had  625  times  fewer  computational  elements,  with 
the  same  20  m  topographic  resolution.  RFSM-EDA  is  capable  of 
further  increasing  the  topographic  resolution,  and  mesh  B  has  the 
same  16  elements,  but  with  a  2m  topographic  resolution,  10 
times  greater  than  the  other  models.  Mesh  C  has  the  recom¬ 
mended  10  000  elements.  See  Table  2  for  details. 

Figure  3  shows  the  model  results  at  test  point  4  (closest  to  the 
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EA  Test  Test  specification 


Details  of  RFSM-EDA  meshes 


Mesh  A 


Mesh  B 


Mesh  C 


Cell  size:  m2 
(~no.  of 
elements) 

Average  IZ 
size:  m2 
(no.  of  IZs) 

Sub-element 
cell  size:m2 

Average  IZ 
size:  m2 
(no.  of  IZs) 

Sub-element 
cell  size:  m2 

Average  IZ 
size:  m2 
(no.  of  IZs) 

Sub-element 
cell  size:  m2 

2A 

400  m2 

250  000  m2 

400  m2 

250  000  m2 

4  m2 

400  m2 

400  m2 

(10  000) 

(16) 

(16) 

(10  000) 

4 

25  m2 

2  300  m2 

25  m2 

Not  undertaken 

Not  undertaken 

(80  000) 

(861) 

5 

2  500  m2 

35  000  m2 

2  500  m2 

Not  undertaken 

2  500  m2 

2  500  m2 

(7  600) 

(530) 

(7  643) 

8A 

4  m2 

212  m2 

4  m2 

174  m2 

0-25  m2 

Not  undertaken 

(97  000) 

(1  786) 

(2  207) 

Table  2.  Details  of  mesh  sizes  and  sub-element  cell  resolutions 
for  meshes  A,  B  &  C. 


Test  2 A  -  point  4  Test  2 A  -  point  5 


Figure  3.  Level  plots  for  test  2  -  points  4  and  5 


boundary  condition)  and  point  5  (the  farthest  point  from  the 
boundary  condition,  which  receives  a  significant  flow  of  water). 
At  point  4,  RFSM-EDA’s  mesh  A  results  have  a  similar  profile  to 
the  other  models,  but  the  peak  level  is  ~4-6  cm  lower.  The  final 
level  matches  the  other  models  exactly.  At  point  5  there  is  a  large 
spread  in  the  results  of  all  models,  not  just  the  ones  shown  here 
(Environment  Agency,  2010).  Even  so,  RFSM-EDA’s  results 
closely  match  those  of  Tuflow  and  Info  Works,  and  from  6  h  on 
they  remain  within  2  cm  of  Tuflow.  RFSM-EDA  predicts  the 
water  levels  to  rise  ~2-5  h  earlier  than  the  other  models.  This  is 


due  to  the  large  IZs  of  RFSM-EDA.  When  water  over-tops  the 
preceding  crest  it  immediately  fills  up  from  the  IZ  base  (the 
location  of  the  test  point).  For  the  other  models  the  water  must 
travel  through  a  number  of  cells  after  the  crest  before  it  reaches 
the  test  location. 

Mesh  B  produces  a  response  that  is  quite  different  from  the  other 
models  (Figure  3).  Although  it  matches  the  peak  level  of  the 
other  models  at  point  4,  the  final  level  is  1-2  cm  higher.  This  is 
because,  rather  than  averaging  the  DEM  to  20  m,  it  utilises  all 


6 


Water  Management 

Volume  165  Issue  WM1 


A  highly  efficient  2D  flood  model  with 
sub-element  topography 

Jamieson,  Lhomme,  Wright  and  Gouldby 


the  topographic  information  available  at  a  2  m  resolution.  There¬ 
fore,  this  mesh  depicts  the  crests  with  a  higher  level  of  accuracy 
than  the  models  using  the  averaged  20  m  DEM.  Once  the  water 
has  spread  over  several  depressions  and  reached  point  5,  there  is 
a  noticeable  cumulative  effect;  the  water  levels  rise  significantly 
only  after  35  h.  For  this  test,  therefore,  the  topography  has  a 
greater  impact  on  the  results  than  the  physical  complexity  of  the 
model. 

RFSM-EDA  is  also  used  with  mesh  C,  which  matches  the  test 
specification  with  10  000  IZs.  As  would  be  expected,  the  results 
have  a  close  match  to  the  other  models.  At  point  4  they  remain 
with  6  mm  of  JFLOW,  and  at  point  5  the  results  lie  in  the  middle 
of  all  the  others,  and  are  closest  to  JFLOW.  Although  not  a  model 
validation,  this  demonstrates  that  RFSM-EDA  behaves  as  ex¬ 
pected  when  used  with  the  same  computational  resolution  as  the 
other  models. 

The  mesh  A  results  show  a  significant  improvement  over  the 
older  dynamic  RFSM,  which  reaches  a  peak  level  approximately 
10  cm  lower  than  the  other  models  at  point  4.  At  point  5  the 
dynamic  RFSM’s  levels  rise  much  too  fast  and  finish  ~4  cm 
higher  than  with  Info  Works  and  Tuflow. 

The  RFSM-EDA  simulations  using  meshes  A  and  B  were 
computationally  fast,  with  equal  runtimes  of  ~0-9  s.  A  large 
proportion  of  this  time  was  spent  communicating  with  the  SQL 
database  that  holds  the  data,  and  therefore  increasing  or  decreas¬ 
ing  the  alpha  value  had  little  or  no  effect  on  simulation  runtimes, 
and  the  increased  topographic  resolution  of  mesh  ‘B’  did  not  slow 
the  model  relative  to  mesh  ‘A’.  The  depression  shape  of  the  IZs 
meant  that  there  was  a  natural  resistance  to  mass  balance  errors. 
The  simulations  were  completed  with  alpha  values  of  1,  with 
median  time  steps  of  ~62  s.  No  instabilities  were  found  and  the 
mass  balance  errors  were  0%.  The  simulation  with  mesh  C  had 
only  one  sub-element  cell  in  each  IZ,  so  did  not  benefit  from  the 
IZ  depression  shape.  However,  it  was  also  able  to  use  an  alpha 
value  of  1  with  only  a  0-3%  mass  balance  error. 

Overall  the  results  of  test  2  show  that  RFSM-EDA  can  effectively 
predict  propagation  of  flood  waters  over  a  complex  domain.  This  is 
encouraging  given  that  only  16  computational  elements  are  used. 

3.3  Test  4:  Rate  of  propagation  over  extended 
floodplains 

The  speed  of  propagation  of  a  flood  wave  is  tested  in  test  4.  A 
completely  flat  domain  is  used,  with  an  inflow  hydrograph  applied 
at  the  centre  of  the  left  boundary,  to  produce  a  semi-circular  flood 
wave.  It  is  not  possible  to  automatically  generate  the  IZs  as  there 
is  no  topographic  variation  in  the  domain;  a  regular  grid  has 
therefore  been  used.  The  specified  resolution  is  5  m  with 
approximately  80  000  elements.  For  this  test,  RFSM-EDA  uses 
mesh  A  with  861  elements,  93  times  fewer  than  the  specification. 
There  is  no  value  in  assessing  results  of  mesh  B  or  C  due  to  the 
flat  topography. 


Figure  4  shows  15  cm  depth  contours  at  1  h  and  3  h  after  start  of 
inundation.  The  coarse  resolution  means  RFSM-EDA  is  not  able 
to  resolve  the  wetting  front  to  the  same  level  of  detail  as  the  other 
models.  However,  the  speed  of  propagation  is  a  significant 
improvement  over  the  dynamic  RFSM,  which  appears  too  slow 
and  also  appears  to  exhibit  some  oscillatory  behaviour  at  the 
wetting  front.  InfoWorks  and  Tuflow  predict  the  flow  boundary  in 
concentric  circles,  whereas  RFSM-EDA  exhibits  a  very  slight 
preferential  flow  towards  the  diagonals,  similar  to  Lisflood-ACC. 
This  has  been  seen  in  several  models  that  have  the  x  and  y  flow 
directions  decoupled  (Neal  et  al.,  2011).  For  all  tests  completed 
by  RFSM-EDA,  this  pattern  has  only  been  observed  on  perfectly 
flat  topography  when  using  a  regular  grid.  It  is  therefore  not 
expected  in  real  topographic  environments. 

Figure  5  shows  the  level  plot  at  point  2,  100  m  from  the  inflow. 
As  the  RFSM-EDA’s  grid  cells  are  large,  the  water  reaches  test 
point  2  marginally  before  the  other  models,  but  the  form  of  the 
curve  matches  those  of  the  other  models  well,  with  a  peak  level 
of  261cm  compared  with  26-6  cm  for  Lisflood-ACC  and 
27-5  cm  for  Tuflow.  RFSM-EDA’s  velocity  profile  has  the  correct 
shape,  although  the  results  are  too  low,  with  a  peak  velocity  of 
0-20  m/s  compared  with  0-23  m/s  and  0-25  m/s  for  Lisflood-ACC 
and  Tuflow,  respectively.  This  may  be  due  to  the  assumptions 
used  in  the  area-averaging  for  the  velocity  calculations. 

RFSM-EDA  completed  the  test  using  an  alpha  value  of  3,  which 
produced  a  runtime  of  ~13  s,  significantly  faster  than  any  other 
model  (the  fastest  other  model  took  nearly  six  times  longer).  This 
was  achieved  with  zero  mass  balance  errors. 

3.4  Test  5:  Valley  flooding 

Test  5  simulates  a  major  flood  inundation  from  a  dam  failure  in  a 
valley.  The  test  domain  has  a  constant  downward  slope  with  a 
hydrograph  applied  at  the  top  of  the  valley.  For  mesh  A,  a  regular 
square  grid  was  adopted  as  few  depressions  could  be  found.  The 
mesh  had  530  elements,  each  200  m  square  (except  at  the  domain 
boundary),  whereas  the  specified  resolution  was  over  14  times 
this  number  of  elements,  with  approximately  7600.  For  mesh  C, 
7643  IZs  where  used  with  one  sub-element  cell  per  IZ  to  match 
that  of  the  other  models. 

Figures  6  and  7  respectively  show  the  first  and  last  test  points  in 
the  domain.  The  IZs  in  mesh  A  generally  have  16  sub-element 
cells  in  them.  The  test  point  will  be  at  one  of  these  cells,  but 
normally  at  least  one  of  the  other  15  sub-element  cells  will  have 
a  lower  level.  This  is  why  the  levels  for  mesh  A  can  be  seen  to 
start  from  a  lower  level  than  the  other  models.  At  point  1  the 
levels  finish  ~28  cm  lower  than  Lisflood-ACC,  and  peak  ~21  cm 
lower.  It  is  clear  that  the  Dynamic  RFSM  did  not  perform  well 
for  these  tests,  which  is  probably  due  to  the  use  of  flow  limiters. 

Mesh  C  shows  that  when  running  at  the  recommended  resolution, 
RFSM-EDA  produces  results  that  are  very  similar  to  the  other 
models.  In  fact,  they  are  almost  indistinguishable  from  those  of 
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Figure  4.  Depth  contours  (15  cm)  at  1  h  (inner  concentric  lines) 
and  3  h  (outer  concentric  lines)  for  test  4 


Lisflood-ACC.  This  is  to  be  expected  as  when  there  is  only  one 
sub-element  cell  in  the  IZ  the  governing  equations  simplify  to  an 
equivalent  of  Bates  et  al.  (2010).  The  differences  seen  between 
RFSM-EDA’s  results  for  mesh  A  and  the  other  models  are 
therefore  primarily  caused  by  the  size  of  the  IZs.  Normally  IZs 
have  a  natural  resistance  to  over-rapid  spreading  of  water,  as  each 
IZ  must  fill  up  a  depression  to  the  crest  level  before  it  can 
continue  to  flux.  However,  on  test  5  there  is  an  almost  constant 


slope  and  very  few  depressions  can  be  found.  The  IZs  fill  up 
from  the  lowest  sub-element  cell,  and  can  immediately  continue 
to  flux,  causing  over-rapid  down-slope  wetting.  This  has  a 
cumulative  effect  down  the  whole  valley.  At  point  1  (~3-2  km 
from  the  inflow)  the  levels  start  to  rise  ~7  min  earlier  than  in  the 
other  models,  but  by  point  5  (~  15-7  km  from  the  inflow)  it  is 
~50  min  too  early.  The  wetting  front  propagates  at  7-1  km/h  for 
RFSM-EDA,  and  at  an  average  5-2  km/h  for  the  other  models. 
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Test  4  -  point  2 


Test  4  -  point  2 


Figure  5.  Levels  and  velocities  for  test  4  -  point  2,  100  m  from 
the  inflow 


The  velocities  predicted  by  RFSM-EDA  (mesh  A)  match  the 
other  models  well.  At  point  1  the  velocities  for  RFSM-EDA 
remain  within  0T5  m/s  of  the  other  models,  except  for  a  ~5  min 
window  at  0-5  h  where  it  peaks  ~0-4- 1  m/s  lower.  Adjusting  for 
the  time  lag,  at  point  5  the  velocities  of  RFSM-EDA  remain 
within  01  m/s  of  the  other  models,  except  for  a  15  min  window 
when  they  are  ~0- 1-0-4  m/s  lower. 

RFSM-EDA  (mesh  A)  was  run  with  an  alpha  value  of  2,  which 
resulted  in  a  final  mass  balance  error  of  only  0  02%.  It  completes 
the  simulation  in  ~14s,  which  is  significantly  faster  than  all 
other  models  that  undertook  the  test  (ranging  from  0-6  to 


Test  5  -  point  1 


Test  5  -  point  1 


Figure  6.  Levels  and  velocities  for  test  5  -  point  1,  near  the 
inflow  source 


350  min).  The  mesh  C  model  was  also  run  with  an  alpha  value  of 
2,  and  had  zero  final  mass  balance  errors. 

3.5  Test  8A:  Rainfall  and  point  source  surface  flow 
This  is  a  test  of  high-resolution  modelling  in  an  urban  environ¬ 
ment,  initially  from  a  global  pluvial  event,  and  subsequently  from 
a  surcharging  culvert  in  the  top  right  comer  of  the  domain.  The 
simulation  is  run  long  enough  to  allow  the  water  to  settle  in  the 
lower  areas. 

This  test  case  has  real  topography  and  RFSM-EDA  can  therefore 
use  its  automatic  mesh  generation.  The  resulting  IZs  have  quite 


9 


Water  Management 

Volume  165  Issue  WM1 


A  highly  efficient  2D  flood  model  with 
sub-element  topography 

Jamieson,  Lhomme,  Wright  and  Gouldby 


Test  5  -  point  5 


Figure  7.  Levels  and  velocities  for  test  5  -  point  5,  at  the  bottom 
of  the  valley 


complex  shapes  and  neighbour  relations,  as  shown  in  Figure  8. 
Mesh  A  has  1786  IZs  with  the  recommended  topographic 
resolution  of  2  m.  This  is  ~54  times  fewer  than  the  specified 
97  000  elements.  Mesh  B  is  also  used,  which  has  2207  elements 
and  a  sub-element  resolution  of  0-5  m. 

The  results  for  RFSM-EDA  are  good  considering  that  the  scale  of 
the  test  is  far  smaller  than  the  model  was  designed  for.  Results 
for  meshes  A  and  B  both  have  levels  approximately  10  cm  higher 
than  the  other  models  at  point  7  for  both  the  first  and  second 
peaks  (Figure  9).  For  point  8  (Figure  10),  mesh  A  results  are 
~8  cm  higher  at  the  first  peak  and  1-2  cm  lower  for  the  second 
peak,  whereas  mesh  B  results  are  1  -2  cm  higher  for  both  peaks. 
The  final  levels  are  similar  for  all  models  in  point  7,  but  very 


widely  spread  for  point  8.  This  indicates  that  the  different  models 
have  likely  sampled  or  averaged  the  raw  DTM  in  different  ways. 
The  dynamic  RFSM  performs  poorly  and  does  not  match  the 
shape  of  the  curves  as  well  as  RFSM-EDA  does.  Although  the 
timing  of  the  velocities  is  good,  the  magnitudes  are  lower  for  both 
meshes;  roughly  half  that  of  the  other  models,  with  mesh  B 
tending  to  have  greater  velocities.  The  IZs  have  complex  shapes 
which  encompass  the  major  flow  routes  on  the  roads  (where  the 
test  points  are  located)  as  well  as  the  areas  surrounding  the  roads. 
It  is  likely  therefore  that  the  lower  velocities  are  a  result  of  the 
velocity  area-averaging  over  the  large  IZs. 

While  the  local  response  of  RFSM-EDA  may  differ  in  a  few 
places  from  the  other  models,  the  overall  model  response  is 
similar,  and  is  illustrated  by  comparing  the  depth  contours  over 
the  domain  (Figure  11);  note  that  to  avoid  complication,  only  the 
results  of  RFSM-EDA  (mesh  A)  and  Tuflow  FV  are  shown  in 
Figure  11.  RFSM-EDA  appears  to  match  Tuflow  very  well  for 
depths  of  20  cm,  but  the  lower  depths  of  5  cm  are  not  very  well 
depicted  in  certain  parts  of  the  domain,  particularly  the  sloped 
areas  to  the  east.  This  is  investigated  further  by  calculating  the 
F-statistic,  which  measures  the  predictive  accuracy  of  the  inun¬ 
dated  area  (Horritt  and  Bates,  2001)  relative  to  the  Tuflow  results. 
Mesh  A  has  an  F  value  of  54%  for  depths  greater  than  5  cm. 
When  the  depth  threshold  is  increased  to  20  cm,  mesh  A  has  an 
F  value  of  69%.  For  mesh  B  the  predictions  are  53%  and  71%, 
respectively.  Clearly  RFSM-EDA  has  some  difficulties  simulating 
the  shallow  flow  paths,  but  when  greater  depths  are  considered  it 
performs  much  better.  There  is  no  major  difference  in  the  results 
of  mesh  A  and  B. 

RFSM-EDA  is  run  with  an  alpha  parameter  of  4  for  these  tests. 
For  mesh  A  this  gives  a  runtime  of  2-90  min,  much  faster  than 
any  other  model,  and  with  a  mass  balance  error  of  only  0'06%. 
For  mesh  B  it  runs  in  4-32  min  with  a  mass  balance  error  of 
0-83%.  The  longer  runtime  for  mesh  B  is  partially  because  it  has 
~24%  more  IZs,  but  it  is  also  due  to  the  finer  sub-element 
resolution.  On  average  an  IZ  in  mesh  B  has  ~27  sub-element 
cells  in  each  interface,  whereas  mesh  A  has  only  ~6.  This  means 
that  simulations  using  mesh  B  have  a  lot  more  calculations  to 
undertake  for  the  interface  fluxes  than  simulations  using  mesh  A. 

3.6  Computational  efficiency 

Unlike  most  similar  models,  the  data  needed  to  run  RFSM-EDA 
is  stored  in  an  SQL  database.  This  allows  for  efficient  modular¬ 
isation  within  probabilistic  modelling  frameworks  like  MDSF2, 
but  it  can  slow  down  the  simulation  through  read  and  write 
access  to  the  SQL  server.  Whilst  this  is  not  generally  an  issue 
unless  a  high  frequency  of  intermediate  results  are  required,  it 
can  dominate  the  performance  in  very  short  tests:  for  example,  in 
test  2A  one-fifth  of  the  simulation  time  is  spent  in  initialisation. 

The  simulation  times  for  all  models/tests  are  given  in  Table  3.  It 
is  important  to  note  that  these  results  may  not  present  a  fair 
comparison,  as  computers  with  varying  specifications  have  been 


10 


Water  Management 

Volume  165  Issue  WM1 


A  highly  efficient  2D  flood  model  with 
sub-element  topography 

Jamieson,  Lhomme,  Wright  and  Gouldby 


Figure  8.  Impact  zones  in  RFSM-EDA  mesh  A  for  test  8A,  with 
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Test  8A  -  point  7 


Test  8A  -  point  7 


Figure  9.  Levels  and  velocities  for  test  8a  -  point  7 


used  and  some  of  the  models  also  used  parallel  processing  (e.g. 
Tuflow-FV,  InfoWorks-ICM  and  JFLOW-GPU). 

Table  3  clearly  shows  that  RFSM-EDA  is  fast;  the  fastest  in  every 
test  attempted.  This  has  been  achieved  without  the  benefit  of 
parallelisation.  It  performs  well  in  these  tests  primarily  because  it 
was  possible  to  maximise  the  benefit  of  the  sub-element  re¬ 
presentation  while  undertaking  the  computations  on  a  coarse  grid. 
At  larger  spatial  scales,  for  which  the  model  has  been  developed, 
further  benefits  are  likely  to  be  realised. 

It  is  also  expected  that  using  a  single  flux  calculation  (based  on 


total  interface  properties  as  opposed  to  the  compound  section 
currently  used)  would  significantly  improve  simulation  runtimes. 

4.  Discussion 

RFSM-EDA  was  designed  to  be  used  on  large  (city/regional) 
scales  with  variable  (i.e.  real)  topography.  The  EA’s  benchmark 
tests  are  small  scale  and  a  number  have  artificially  smooth 
topography.  Despite  this,  the  RFSM-EDA  has  demonstrated  an 
ability  to  generate  results  that  are  in  line  with  those  of  models 
that  comprise  a  more  complex  representation  of  the  physical 
processes  and  thus  take  a  longer  computational  time.  Moreover, 
RFSM-EDA  can  incorporate  even  finer-scale  topography  with 
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Test  8A  -  point  8 


Test  8A  -  point  8 


Figure  10.  Levels  and  velocities  for  test  8a  -  point  8 
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Figure  11.  Depth  contours  for  RFSM-EDA  and  Tuflow  FV  for  test 
8A 


minimal  impact  on  runtimes.  This  reduces  the  need  of  the 
modeller  to  introduce  additional  uncertainties  to  the  modelling 
process  by  averaging  or  re-sampling  the  DEM.  Given  that  much 
of  the  flood  modelling  undertaken  in  the  UK  is  at  larger  spatial 
scales  and  often  of  a  probabilistic  nature,  it  may  be  appropriate  to 
consider  the  introduction  of  additional  tests  that  are  able  to 
appropriately  verify  models  that  are  developed  for  this  purpose. 

The  schematisation  of  RFSM-EDA  means  that  water  fills  from 
the  lowest  point  in  an  impact  zone.  On  the  relatively  rare 
occasion  that  natural  floodplain  depressions  do  not  exist,  such  as 
in  test  5,  water  leaves  an  impact  zone  immediately  upon  wetting. 
This  results  in  an  overestimation  of  propagation  speed  by  36%. 


For  large-scale  probabilistic  modelling  this  source  of  error  is 
unlikely  to  be  significant,  and  the  results  presented  herein  show 
that  peak  levels  and  flows  are  predicted  reasonably  accurately 
(see  Figures  5  and  6).  In  other  situations  where  RFSM-EDA 
results  differ  more  markedly  from  those  of  the  other  models 
presented  here,  it  is  worth  considering  that  the  EA  benchmarking 
report  includes  many  more  model  results  (Wright  et  al,  2012), 
and  the  peak  flood  levels  of  RFSM-EDA  are  within  the  spread  of 
these  results. 

The  sub-element  representation  in  RFSM-EDA  offers  an  effective 
approach  for  reducing  runtime  while  preserving  or,  in  some 
cases,  increasing  topographic  accuracy.  The  results  using  the 
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Model  Computation  time  in  minutes  for  each  test 


2 

4 

5 

8A 

RFSM-EDA  (mesh  A) 

0015 

0-21 

0-23 

2-9 

Dynamic  RFSM 

0-19 

5-8 

9-8 

23-3 

Tuflow  FV 

2-64 

24-5 

2-9 

72-6 

InfoWorks  ICMa 

0-73 

6-5 

0-7 

27  1 

JFLOW-GPU 

1  83 

2-3 

1 0-2 

1 6  2 

Lisflood-ACC  b 

n/a 

1  97 

0-68 

n/a 

Fastest  otherc 

0-4 

1  27 

0-6 

4 

Slowest  otherc 

130 

282-8 

350 

307-8 

a  The  runtimes  are  taken  from  InfoWorks  RS,  but  the  results  in  this  paper  are  from 
InfoWorks  ICM;  little  difference  is  expected. 
b  Lisflood-ACC  runtimes  appear  in  Neal  et  al.  (201 1). 

c  Fastest  and  slowest  models  other  than  those  shown  in  this  paper,  but  appearing  in 
Wright  et  al.  (2012). 


Table  3.  Simulation  runtimes  for  different  models,  fastest  in  bold 
type 


recommended  DEM  resolution  matched  the  other  models  well.  In 
some  cases,  such  as  test  2A,  using  an  even  higher  DEM 
resolution  produced  a  step  change  in  model  response,  which 
implies  that  the  topography  has  a  greater  effect  on  simulation 
results  than  process  representation.  Additionally,  because  the 
mesh  is  automatically  aligned  to  topographic  features  such  as 
embankments  and  dykes,  it  will  always  respect  the  effect  they 
have  on  propagation  directions,  regardless  of  grid  scale. 

The  adaptive  time  step  used  by  RFSM-EDA  has  shown  to  be 
effective  for  all  tests.  Unlike  Lisflood-ACC,  which  generally 
needs  alpha  values  significantly  below  unity  (Neal  et  al.,  2011), 
RFSM-EDA  is  stable  with  a  value  of  1  or  significantly  above. 
The  fact  that  the  alpha  value  could  be  as  large  as  4  in  test  8A 
implies  that  the  CFL  condition  used  (Equation  7)  may  be 
conservative  for  this  algorithm.  It  is  likely  that  this  is  due  to  the 
inclusion  of  the  velocity  vector  in  the  CFL  condition,  which  is 
not  included  in  the  original  model  of  Bates  et  al.  (2010). 
Although  several  other  diffusive  models  use  velocity  in  their 
stability  condition,  such  as  Bradbrook  et  al.  (2004),  an  alternative 
formulation  that  excludes  velocity  may  be  more  appropriate  for 
RFSM-EDA. 

Although  the  results  have  already  been  shown  to  be  good  for 
these  small  tests,  there  is  potential  for  further  improvements.  In 
test  5  the  propagation  speeds  are  too  fast  down  the  valley,  which 
is  primarily  caused  by  the  large  computational  elements.  Future 
work  should  aim  to  find  an  approach  to  limit  the  propagation 
speeds  for  large  computational  elements.  Using  a  single  flux 
calculation  at  the  interface,  rather  than  a  summation  of  panel 
fluxes,  has  the  potential  to  make  RFSM-EDA  considerably  faster 
still,  although  the  impact  on  simulation  accuracy  will  require 


verification.  Some  investigations  may  be  necessary  to  see  whether 
predictions  of  low-depth  flow  paths  can  be  improved,  as  in  test 
8A,  but  these  shallow  flow  paths  are  less  important  for  probabil¬ 
istic  risk  calculations  than  greater  depths. 

The  RFSM-EDA  has  been  developed  specifically  for  use  at  larger 
spatial  scales  and  within  the  context  of  probabilistic  simulations. 
The  model  provides  a  step-change  in  accuracy  over  previous 
versions,  the  dynamic  RFSM  and  the  direct  RFSM  (which  is 
currently  used  within  the  Environment  Agency’s  NaFRA  and 
MDSF2  systems).  This  significant  improvement  comes  with  the 
price  of  additional  computational  expense  over  the  direct  RFSM. 
The  computational  expense  is  however,  a  fraction  of  that 
associated  with  alternative  models  that  solve  the  full  SWE  on 
conventional  grid  systems.  The  model  therefore  provides  a  good 
compromise  between  practical  computational  times,  while  provid¬ 
ing  robust  flood  simulations. 

5.  Conclusions 

RFSM-EDA  has  been  applied  to  six  of  the  EA’s  hydraulic 
benchmarking  tests,  four  of  which  are  shown  in  this  paper.  The 
model  was  designed  to  be  used  on  larger  domains  of  naturally 
varying  topography,  but  nonetheless  has  performed  well  given  the 
small-scale  nature  of  the  tests.  The  peak  levels  predicted  by 
RFSM-EDA  differed  by  less  than  ±10  cm  from  the  other  models 
in  all  cases  except  for  test  5,  where  they  were  within  ±50  cm. 
This  is  a  mesh  effect  rather  than  a  numerical  inaccuracy,  as  when 
using  an  equivalent  mesh  resolution  the  results  were  visually 
identical  to  those  of  Lisflood-ACC.  The  velocity  predictions  had 
a  similar  form  to  the  other  models,  though  they  tended  to  be  20- 
60%  lower.  This  is  primarily  because  an  impact  zone  average 
velocity  is  used;  using  a  maximum  velocity  would  be  more 
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conservative.  RFSM-EDA  clearly  offers  a  step  change  in  accu¬ 
racy  over  the  direct  RFSM  and  dynamic  RFSM,  while  comparing 
favourably  with  industry  standard  codes.  As  RFSM-EDA  can 
increase  topographic  resolution  without  needing  to  increase  the 
number  of  computational  elements,  it  is  able  to  improve  simula¬ 
tion  accuracy  further  with  minimal  change  in  computational 
burden. 

RFSM-EDA  was  the  fastest  of  all  models  by  a  considerable 
margin  on  all  of  the  tests  (less  than  a  tenth  of  the  average  runtime 
of  the  models  shown  here,  and  between  4%  and  73%,  depending 
on  the  test,  of  the  runtime  of  the  otherwise  fastest  model).  It  has 
the  potential  to  be  even  faster  if  simpler  flux  calculations  and 
parallelisation  are  implemented.  Additional  testing  on  very  large 
regional  domains  is  underway,  and  it  is  likely  that  the  benefits  of 
the  scheme  will  become  even  more  apparent  as  the  trade-off 
between  simulation  time  and  grid  resolution  becomes  more  severe 
for  conventional  models. 


RFSM-EDA  has  completed  a  selection  of  the  EA  benchmarking 
tests  with  fast  runtimes  and  results  accurate  enough  for  broad- 
scale  flood  risk  assessments.  The  tests  present  a  proof  of  concept, 
and  demonstrate  that  the  model  has  the  potential  to  be  an 
effective  tool  for  large-scale  and  probabilistic  inundation  model¬ 
ling. 
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